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ABSTRACT 

In this work, the feasibility of using massively parallel computation to study the 
response of ablative materials is investigated. Explicit and implicit finite difference methods 
are used on a massively parallel computer, the Thinking Machines CM-5. The governing 
equations are a set of nonlinear partial differential equations. The governing equations are 
developed for three sample problems: (1) transpiration cooling, (2) ablative composite plate, 
and (3) restrained thermal growth testing. The transpiration cooling problem is solved using 
a solution scheme based solely on the explicit finite difference method. The results are 
compared with available analytical steady-state through-thickness temperature and pressure 
distributions and good agreement between the numerical and analytical solutions is found. 
It is also found that a solution scheme based on the explicit finite difference method has the 
following advantages: incorporates complex physics easily, results in a simple algorithm, and 
is easily parallelizable. However, a solution scheme of this kind needs very small time steps 
to maintain stability. A solution scheme based on the implicit finite difference method has 
the advantage that it does not require very small times steps to maintain stability. However, 
this kind of solution scheme has the disadvantages that complex physics cannot be easily 
incorporated into the algorithm and that the solution scheme is difficult to parallelize. A 
hybrid solution scheme is then developed to combine the strengths of the explicit and implicit 
finite difference methods and minimize their weaknesses. This is achieved by identifying the 
critical time scale associated with the governing equations and applying the appropriate 
finite difference method according to this critical time scale. The hybrid solution scheme is 
then applied to the ablative composite plate and restrained thermal growth problems. The 
gas storage term is included in the explicit pressure calculation of both problems. Results 
from ablative composite plate problems are compared with previous numerical results which 
did not include the gas storage term. It is found that the through-thickness temperature 
distribution is not affected much by the gas storage term. However, the through-thickness 
pressure and stress distributions, and the extent of chemical reactions are different from the 
previous numerical results. Two types of chemical reaction models are used in the restrained 
thermal growth testing problem: (1) pressure-independent Arrhenius type rate equations and 
(2) pressure-dependent Arrhenius type rate equations. The numerical results are compared 
to experimental results and the pressure-dependent model is able to capture the trend better 
than the pressure-independent one. Finally, a performance study is done on the hybrid 
algorithm using the ablative composite plate problem. It is found that there is a good 
speedup of performance on the CM-5. For 32 CPUs, the speedup of performance is 20. The 
efficiency of the algorithm is found to be a function of the size and execution time of a given 
problem and the effective parallelization of the algorithm. It also seems that there is an 
optimum number of CPUs to use for a given problem. 
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CHAPTER 1 


INTRODUCTION 


Ablative composite materials have been used in a broad range of 
applications. The production of solid fuel rocket motors, planetary-entry 
probes, and reentry vehicles was made possible by these materials. Ablatives 
are used to protect structures from extreme temperature environments. 
When these materials are exposed to a high heat flux environment, extreme 
thermal gradients, internal pressures due to chemical reactions, and thermal 
and mechanical stresses all develop, which can cause premature failure. 
Hence, to make full use of these materials, all of these aspects of their 
response must be more completely understood. 

The physical phenomena that happen inside these materials can be 
characterized by the following processes. When the surface of the ablative is 
exposed to a high temperature environment, heat is conducted into the 
material. The temperature of the material below the surface will then rise 
and a temperature gradient is established inside the material. When the 
temperature inside the material reaches a high enough value, chemical 
reactions take place. Gases are generated due to these reactions. Due to the 
relatively low permeability of the material, these gases are trapped in the 
material and cause internal pressures to develop. Both the temperature 
gradient and internal pressure will cause stresses to develop in the materials. 
Sometimes the values of these stresses are high enough to cause premature 
failure of the ablative materials. The response of the ablative materials 
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undergoing these processes needs to be understood to prevent premature 
failure. 

There are typically two approaches to study the responses of composite 
ablatives. The experimental approach usually requires large scale test 
specimens which can be prohibitively expensive to manufacture. Moreover, 
these experiments do not always reveal the details of the underlying physical 
processes. The analytical approach can give insight to the physical processes 
that lead to premature failure of such materials if the modeling is done 
properly. However, in order to obtain these insights, the nonlinear partial 
differential equations of a very complete analytical model need to be solved. 
Closed-form solutions are impossible to obtain in most cases. Therefore, 
numerical solutions are needed. 

Two major numerical solution schemes that have been applied to this 
problem are the finite element method (FEM) and the finite difference 
method (FDM). In both methods, some simplifications must be made in the 
governing equations in order to keep the computation tractable. However, it 
is important that sufficient complexities are included in the computations so 
that the predicted results can be used with confidence. A great deal of 
computational power is needed to allow such complexities to be included in a 
solution algorithm. With the recent arrival of massively parallel computers 
such as the Thinking Machines' Connection Machine 5 (CM-5), enough 
computational power has become available to both greatly increase the 
accuracy of such solutions and lower their turn-around time. 

However, to make efficient uses of massively parallel computers, an 
appropriate solution algorithm needs to be developed. Relatively speaking, 
FDM results in simpler algorithms than FEM on a parallel computer. 
Simpler algorithms allow more complexities to be incorporated into an FDM 
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algorithm. Results computed based on this algorithm will simulate reality 
more closely. Therefore, a solution algorithm based on FDM is developed on 
the CM-5 in this work. 

In FDM, there are two major types of solutions scheme. They are the 
explicit finite difference method (EFDM) and the implicit finite difference 
method (IFDM). The two schemes differ in the way they approximate the 
derivatives in the governing equations. There are two major advantages of 
using EFDM. The first advantage is that it leads to a simple algorithm which 
is easy to program. This implies that complex nonlinear physics is relatively 
easy to incorporate into the program. The second advantage is that EFDM is 
well suited to parallel computation. At each time step, the difference 
equations can be solved at all nodes simultaneously and the solution of the 
equations at each node requires knowledge only of the states of its immediate 
neighboring nodes. This minimizes the required communication between 
processors, keeping the parallel computation efficient. 

The major drawback of the EFDM is that computations need to be done 
at relatively small time steps to maintain stability. The allowable time steps 
are usually limited by the fastest physical process associated with the 
problem. In the ablative case, the fastest physical processes are the 
deformations of the solid material and the internal flow of gases. These two 
processes can limit the time step to as small as 10 -6 second. Therefore, 
modeling using a fully-explicit scheme may be impractical for a typical 
simulation time which is on the order of 100 seconds. A hybrid algorithm 
which will be discussed later is proposed to alleviate the time-step problem. 

In this work, approaches to solve ablation type problems on massively 
parallel computers are explored. Analytical models and numerical solution 
schemes are developed to solve several physical problems of interest. 
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The three physical problems considered are: (1) transpiration cooling of 
a plate, (2) ablation of a composite plate, and (3) restrained thermal growth 
testing. The transpiration cooling problem is used to verify the EFDM. The 
internal pressure, temperature, and stress distributions are obtained using 
EFDM and compared to available analytical solutions. The ablative 
composite plate problem is used to demonstrate the capability of the EFDM 
scheme to incorporate complex physics and solve parctical problems. A 
typical solid rocket motor nozzle liner problem is solved, and the solution 
compared to previous numerical solutions of the problem. The restrained 
thermal growth problem is used to demonstrate the adaptability of the 
numerical method. Parametric studies are performed on the effects of 
different chemical reaction models and the results compared with 
experimental data. 

The present work is organized as follows. In Chapter 2, previous work 
and relevant background on massively parallel computing are discussed. A 
concise problem statement is given in Chapter 3. In Chapter 4, the general 
governing equations are developed. The governing equations are then 
simplified appropriately for each of the three cases studied: transpiration 
cooling, ablative composite plate, and restrained thermal growth. The results 
are presented in Chapter 5. The performance of the hybrid algorithm on the 
CM-5 is also presented in Chapter 5. The conclusions drawn from the results 
and recommendations for future research are presented in Chapter 6. 
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CHAPTER 2 


BACKGROUND 


The first section of this chapter contains a discussion of the previous 
analytical and experimental work done in the area of ablative composite 
materials. A brief discussion of previous work on numerical schemes then 
follows. In the second section of this chapter, relevant information 
concerning parallel computation on machines such as the CM-5 is presented. 
The first part of this section contains a brief discussion on the architecture of 
the CM-5, as well as how codes using CM-FORTRAN [1] should be written to 
take advantage of the architecture. In the second part of this section, the 
metrics for measuring the performance of parallelized codes is presented. 

2.1 Experimental Investigations and Modeling 

The earliest work done on natural ablative materials was done on wood 
for fire-retarding purposes (see [2] for more early history). For man-made 
ablatives, one of the earlier works was done by Moyer and Rindal [3] in 1963 
in which the thermal response of materials used as heat shields for reentry 
vehicles was investigated. 

Henderson and his colleagues did extensive experimental work to 
determine the properties of glass-phenolic ablative composite materials in the 
early 1980's [4]. Stokes [5] and Hubbert [61 performed many experiments to 
study the material response of ablative composite materials. The one of 
particular interest to this study is the restrained thermal growth (RTG) test. 
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McManus and Springer [7] also did some experimental work on carbon- 
phenolic ablative materials to validate his model and the CHAR computer 
code. Florio et al. [8] experimentally determined the volumetric heat transfer 
coefficient in decomposing polymer composites. This volumetric heat transfer 
coefficient was used in their study of the assumption of local thermal 
equilibrium [9]. 

Much work has been done on modeling the behavior of ablative 
composite materials. Henderson developed a simple model which included 
an Arrhenius reaction model, an energy equation, and a steady-state mass 
flow equation [4]. The model was later refined to include a mass flow 
equation based on Darcy's law and the thermal expansion of the solid 
material [10]. These models predicted the internal pressure and temperature 
predictions distributions but did not give stress distributions. Kuhlmann 
[11], McManus [12], Sullivan [13], and Weiler [14, 15] developed models 
which also included the stress distributions inside the solid material. This 
was achieved by applying the theory of poroelasticity [16-21] in combination 
with the existing thermochemical and gas flow theories to predict the 
material's temperature, chemical state, internal pressure, and stresses 
simultaneously. Each author used a different approach to derive the coupled 
thermoporoelastic equations. These governing equations are highly 
nonlinear due to the fact that the coefficients in the constitutive relations are 
functions of the independent variables (temperatures, pressures, stresses, 
etc.). When complex chemical reaction models are used, governing equations 
can be made even more nonlinear. For example, Tai [22] developed a new 
evaporation model for ablative composite materials which is a function of 
temperature and pressure. 
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Sullivan [23, 24] proposed a thermodynamic approach to derive a new 
set of constitutive relations for decomposing ablative materials. In general, 
the coefficients in the newly developed constitutive relations are functions of 
the independent variables and thus the governing equations are highly 
nonlinear. This approach was proposed to overcome the limitations of the 
previous models which were all based on porous media flow and 
poroelasticity. These limitations are: (1) gases generated from chemical 
reactions act together as a single equivalent fluid, (2) a well-defined boundary 
exists between the fluid and solid constituents where mechanical equilibrium 
is maintained and (3) the forces that exists between the solid and fluid 
constituents are purely mechanical in nature. However, certain carbon fibers 
used in polymeric composites are known to be hydrophilic (i.e. attract water) 
when heated to high temperatures [251. This is due to the presence of 
activated carbon sites on their surface. Chemical forces then develop between 
the carbon fibers and water molecules liberated from the resin. These forces 
may have significant influence on the mechanical behavior of the material. 
In the only work to date on this problem, chemisorption of H 2 0 in the matrix 
was dealt with on an ad-hoc basis in the model developed by Tai and 
McManus [22]. 

2.2 Numerical Methods 

Due to the complex physics occurring inside an ablative composite 
material, the models developed so far require the solution of a set of highly 
nonlinear partial differential equations. Closed form solutions for these 
governing equations are impossible to obtain. Therefore, numerical schemes 
have been used. Two types of numerical schemes have been adopted by 
researchers in this field. Sullivan, Kuhlmann, and Weiler adopted the finite 
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element method. Other researchers such as Henderson and McManus have 
used the finite difference method. The finite element method has the 
advantage that geometry and boundary conditions can be accurately modeled. 
The major disadvantage of the finite element method is the difficulty in 
incorporating complex physics (nonlinear constitutive laws for example). The 
advantage of the finite difference method is that it is relatively easy to 
incorporate complex physics. The disadvantage of the finite difference 
method is that geometry and boundary conditions are harder to model than 
in the finite element method. Within the finite difference method, there are 
two ways a differential equation can be approximated: (1) the explicit finite 
difference method (EFDM) and (2) the implicit finite difference method 
(IFDM). It is easier to incorporate complex physics into the EFDM than the 
IFDM. However, the EFDM generally requires more computational power 
than the IFDM. Henderson's original work used the EFDM, because the 
EFDM was relatively easy to implement [4]. However, it was later 
abandoned due to limited computational power available at that time. With 
the recent advances in computational power, the explicit finite difference 
method is again becoming a viable method for solving the governing 
equations. 

In this work, a numerical scheme called the hybrid algorithm is 
developed to solve the governing equations. The hybrid algorithm uses both 
the EFDM and the IFDM to solve the governing equations. For this reason, a 
brief discussion of these methods is presented using an example based on a 
simple 1-D heat transfer equation. 

The 1-D heat transfer equation is: 


cTT d 2 T 

dt dx 2 


( 2 - 1 ) 
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The EFDM, using forward difference in time and central difference in space 
[25], can be used to derive the following difference equation: 

T>* x _ T> T> j. T> 

= k \ (2-2) 

A t Ax' 

where Tj is the temperature at the yth time step and the ith node, At is the 
time step, and Ax is the spacing between two nodes. By rearranging Eqn. 

2-2, an expression for T'*' is obtained: 


7f 1 = rTj_ { + (1- 2 r)T{ + rT^ (2-3) 

where r = (kAt)/Ax 2 . It can be seen from Eqn. (2-3) that the value of 
temperature at the next time step for a given node is found from the known 
current values of temperature of that node and its immediate neighbors. A 
finite difference method where the unknown values can be expressed directly 
in terms of the known values is called an EFDM. The scheme used in Eqn. 
(2-3) is illustrated graphically in Figure 2-1. In Figure 2-1, the y axis 
represents time and x axis represents nodal position. As shown in Figure 2- 
1, the value of temperature at the ith nodal point and (y + l)th time step is 
updated by the known values of the neighboring nodes ((i-l)th, ith, and 
(i + 1 )th nodes) at the jth time step. Eqn. 2-3 is stable for time step sizes that 
are less than (Ax 2 /2k ) [25]. 

To illustrate IFDM, the Crank-Nicholson scheme is used [26]. In this 
scheme, the same difference technique as the one in EFDM is used. The only 
difference is that the spatial derivatives are approximated by taking the 
average of its central difference at the jth. and j + 1th time step. By applying 
the IFDM to Eqn. 2-1 the resulting difference equation is: 


j _ y\ 


At 


j__ _ _ 


kf 77 

2 ,' 


,«,-27y+r'_, tx - ur + t>:/ 

Ax 2 


Ax 2 


(2-4) 
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Unknown Value of T 



Figure 2-1. Explicit finite difference method for the simple heat conduction 
equation (Eqn. 2-3). 
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Rearranging Eqn. 2-4 into a more convenient form, 



-(1 + r)7T 


+ -T£ 
2 1 




(2-5) 


It can be seen readily from Eqn. 2-5 that the value of temperature at a given 
node is dependent on both the known (jth time step) and the unknown 
O' + 1th time step) values of temperature at that and neighboring nodes. The 
value Tf*' cannot be solved directly from Eqn. 2-5 as in Eqn. 2-3. If there are 
N internal nodal points then at the y' + lth time step Eqn. 2-5 gives N 
simultaneous equations for the N unknown values in terms of the boundary 
and y th time step values. Such a difference method is called the IFDM. The 
IFDM used in Eqn. 2-5 is illustrated graphically in Figure 2-2. The value 
T/ +l depends on both the current ( jth time step) and future ( j + 1th time step) 
value of temperature at that and neighboring nodes. The future value of the 
neighboring nodes depend on the values of their neighbors, and so on, until 
the problem becomes fully coupled. 

If the coefficient ( k in our example) is itself a function of temperature, 
Eqn. 2-1 becomes non-linear. This presents no complication to solving the 
EFDM (Eqn. 2-3), as k is calculated at timestep j and is thus known. 
Equation 2-5 requires values of k at both timestep j and timestep j+ 1; thus 
the matrix solution becomes nonlinear and must be solved iteratively. When 
k is a constant, Eqn. 2-5 is stable for all positive values of time step size [26]. 
However, a reasonable time step size must be used to maintain accuracy. 

The key characteristics of EFDM and IFDM are shown in the above 
example. In EFDM, the solution can be obtained from known neighboring 
nodal values. However, the time step size must be smaller than a given value 
for EFDM to remain stable. A matrix solution is needed in EFDM and it may 
need to be iterated to solve non-linear problems. IFDM is more stable than 
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Figure 2-2. Implicit finite difference method (Crank-Nicholson Scheme) for 
the simple heat conduction equation (Eqn. 2-3). 
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EFDM. In the above example, the IFDM is stable for all positive values of 
time step size when k is constant . EFDM is more attractive than IFDM for 
implementation on parallel computers, since the algorithm can be 
parallelized more easily (the solution at each node can be obtained directly 
from known neighboring nodal values). Although it is easier to implement 
EFDM on parallel computer, sometimes the time step size necessary for 
stability may become too small for EFDM to be practical. In that case, one 
may need to use IFDM. 

2.3 Parallel Computing on the CM-5 

The Connection Machines CM-5 is a massively parallel, SIMD (Single- 
Instruction Multiple-Data) and MIMD (Multiple-Instruction Multiple-Data) 
computer. Machines of this type consist of a very large number of processing 
elements. Each parallel processing element has its own physically connected 
memory. Intense communication takes place between the processors when 
data needs to be moved from one memory to another. Efficient algorithms 
will minimize this communication. For large numerical codes, the regularity 
of the data structure and the absence of sequential operations are important 
factors to minimi ze interprocessor communications on the CM-5 [27]. 

A schematic drawing of the CM architecture is shown in Figure 2-3. 
The serial control processor directs the actions of a set of parallel processors. 
The serial controller also performs all sequential operations. During such 
operations, the parallel processors do nothing. The parallel processors act on 
data elements stored in their local memories. Parallel processors are most 
efficient when acting on large data set, each element of which can be acted on 
independently. The entire set can then be acted on simultaneously, one 
processor working on each element. 
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Serial Control Parallel 

Computer Processors 



Local Memories 


Figure 2-3. Connection Machine architecture [1]. 
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CM FORTRAN is the language used here to implement the solution 
algorithm on the CM-5 machine. It allows data parallel programming in a 
language familiar to most researchers since it is actually based on FORTRAN 
90. A full description of CM FORTRAN may be found in the Connection 
Machine documentation [1]. 

CM FORTRAN does not require the programmer to be concerned about 
the details of parallelization. The CM FORTRAN compiler performs the 
parallelization after the code is written. However, the programmer needs to 
arrange the data structure so that the compiler can parallelize the code in the 
most optimal way. Optimal parallelization can be achieved by the compiler 
when data structure is arranged into different sets of conformable arrays 
(arrays that are the same size and shape) and all operations are performed 
with conformable arrays from the same set. This allows the compiler to 
assign conformable arrays to the same parallel processor set. When this is 
done, operations are performed in parallel without communications between 
different sets of processors. 

With the above ideas in mind, it can be seen that the EFDM is a 
suitable method to solve differential equations on the CM-5. Take the simple 
heat conduction equation (Eqn. 2-1) as an example. Eqn. 2-3 is obtained 
using EFDM. For demonstration purpose, it is assumed that there are five 
nodal points in the mesh. Using the EOSHIFT function, a utility in CM 
FORTRAN that allows the location of array elements to be shifted by a 
specified amount, three independent conformable arrays containing the 
values of, 7/, 7/_,, and 7/ +) at all nodes can be formed. The EOSHIFT 
function inserts specified values in the appropriate end of an array which are 
used to define boundary conditions. Then the values of Tj* 1 at ail nodes can 
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be obtained with a single operation (28). This process is illustrated in Figure 
2-4. For this reason, the EFDM is very efficient for parallel computations. 

The opposite is true for the IFDM. As shown in Eqn. 2-5, the 
unknowns values of temperature are coupled together. Therefore, in order to 
obtain solutions, a system of equations needs to be solved at the same time. 
Typical numerical schemes available for solving systems of equations (Gauss 
elimination, LU decomposition, and Gauss-Seidel) involve mainly sequential 
operations which require intense interprocessor communication during which 
the parallel processors remain idle. Therefore, the computation becomes less 
efficient. 

Although the EFDM is a very efficient method for solving differential 
equations on the CM-5, it does have a very stringent stability criterion. As 
mentioned before in chapter 1, the time step for the ablative problem can be 
as small as KT 6 second. The IFDM has a much less stringent stability 
criterion than the EFDM. In the 1-D heat conduction example, the IFDM is 
actually unconditionally stable. However, the implementation of IFDM on 
the CM-5 is not as efficient as that of EFDM. 

For any given problem, it is difficult to predict o priori which method is 
the most suitable. The most efficient algorithm will take advantage of the 
strengths of both methods while minimizing the weaknesses. The hybrid 
algorithm is developed based on this concept. 

2.4 Parallel Computing Performance Measures 

For serial algorithms, performance measurement based on millions of 
floating point instructions per second (MFLOPS) is usually appropriate. 
However, applying this type of measurement to a parallel algorithm is not 
appropriate. The reasons are: (1) extra work is done by a 
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Parallel Processors 



* Values of T' and T ' 6 depend on boundary conditions 


Figure 2-4. 


Example computation using parallel processors. 
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parallel computer in the background and (2) synchronization and 
communication overhead costs are not reflected in FLOPS. Moreover, it does 
not give any measure on how effectively the code is parallelized. 

The performance of parallel algorithms is most commonly measured in 
terms of speedup. Speedup of an algorithm executed using N processors is 
defined as: 

S = — (2-6) 

b 

where 5 is the speedup, r, is the execution time using one processor, and t N is 
the execution time using N processors. 

However, simply measuring the speedup, S, is not sufficient to learn 
how well a parallel code is written. By applying Amdahl's law [29], one can 
measure, ideally, how much of an algorithm cannot be parallelized and must 
be run sequentially. This measure is provided by Amdahl's fraction, a. The 
idealized t N can then be written in terms of a and r, as: 

r„ = ar,+(l + a)-^ (2-7) 

Substituting Eqn. (2-7) into Eqn. (2-6) and letting N approach infinity, the 
maximum ideal speedup is obtained: 

S = - (2-8) 

N-tm QT 

The maximum ideal speed up approaches asymptotically to a number 
governed by a which is the fraction of the sequential algorithm that cannot 
be parallelized. 

It is assumed that the Amdahl's fraction, a, is a constant and depends 
only on the algorithm. In most engineering problems, a depends not only on 
the algorithm but also on the problem size. Hence, Amdahl's law is not 
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directly applicable to measure the performance of an algorithm on CM-5, as 
overhead costs such as initial setup, opening and closing files, input/output of 
results, interprocessor communication, and synchronization delays are not 
taken into account. 

To overcome the limitations in Amdahl's law, a measure called the 
effective parallelization, p, has been introduced by E.J. Plaskacz et. al. [30]. 
In order to compute p, the equivalent Amdahl's fraction, a t , is calculated 
first from the measured speedup vising Eqns. (2-6) and (2-7): 




(2-9) 


Note that a e represents the fraction of the code that is running serially, and 
this includes the serial part of an algorithm as well as the overhead costs. 
The effective parallelization, p, is then given by 


p = 1- a t (2-10) 

where p represents the fraction of the code that can be completely 
parallelized. 

Efficiency, e, and excess time,/„ , are two additional useful measures of 
a parallel code's performance. Efficiency is defined by 


e - 


S_ 

N 


( 2 - 11 ) 


Excess time measures the time spent by a processor over and above the time 
required due to ideal speedup and it is given by the following equation: 


«- 12 ) 

The excess time provides a measure of the total time a parallel computer 
spends on overhead costs. 
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CHAPTER 3 


PROBLEM STATEMENT 


In this work, analytical models and numerical solution schemes are 
developed to solve three ablation-type problems of interest. The three 
problems are: (1) transpiration cooling of a uniform plate (2) ablation of a 
composite plate, and (3) restrained thermal growth of a composite test 
specimen. 

In the transpiration cooling problem, a flat plate made of porous 
material with gas flowing through the thickness is considered. The in-plane 
dimensions of the plate are much greater than the through-thickness 
dimension and all boundary conditions are uniformly applied over the in- 
plane dimensions of the plate. A one-dimensional analysis in the thickness 
direction is developed. Given the material and flowing gas properties 
(assumed constant), the initial conditions (temperature and pressure), and 
the boundary conditions (temperature and pressure values at both 
boundaries, tractions at one end, and displacements at the other), 
temperature, pressure, and stress as functions of time and position through 
the thickness, and mass flux as a function of time, are obtained. The solution 
of this problem is used for comparison to exact steady state solutions, and to 
explore the time scales of different physical aspects of the problem. 

In the ablative composite plate problem, a flat plate made of carbon- 
phenolic material, exposed on one surface to a high temperature environment 
and insulated on the other, is considered. The in-plane dimensions of the 
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plate are much greater than the through-thickness dimension and all 
boundary conditions are uniformly applied over the in-plane dimensions of 
the plate. A one-dimensional analysis in the thickness direction is developed. 
Given the material and gas properties, and the initial and boundary 
conditions, temperature, pressure, and stress as functions of time and 
position through the thickness, and the maximum pressure as a function of 
time, are obtained. Initial conditions specified are temperature and pressure, 
uniformly distributed through the thickness of the plate. Boundary 
conditions specified on the exposed surface of the plate are the heat flux and 
ambient pressure values. On the insulated surface, the heat flux and gas 
mass flux values are set to zero. The solutions of this problem are used to 
demonstrate the capability of the massively parallel computer algorithm to 
incorporate complex physics, explore the effects of including additional 
physics on the solutions, and study the performance of the algorithm. 

In the restrained thermal growth (RTG) problem, a cylindrical 
specimen made of carbon phenolic material heated uniformly at a constant 
rate and held at a constant longitudinal strain is considered. A one- 
dimensional model, approximating the cylindrical geometry of the specimen 
as a strip, is developed. Given the material and gas properties and the initial 
and boundary conditions, the restraining stress required to hold the specimen 
at a constant longitudinal strain is obtained as a function of temperature. 
The properties of the material and gas are functions of temperature and 
chemical state. The initial conditions are uniformly distributed temperature 
and pressure values. The boundary conditions at one end of the strip are 
the values of ambient pressure and surface tractions. At the other end of the 
strip, the values of the mass flux and displacements are set equal to zero. 
The solution of the problem is used to demonstrate the capability of the 
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computer algorithm to incorporate complex physics, to perform a parametric 
study of two chemical reaction models, and to compare the analytical results 
to experimental measurements. 
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CHAPTER 4 


THEORY AND IMPLEMENTATION 


In this chapter, the general governing equations are developed for 
three physical problems. Explicit finite difference method (EFDM) or hybrid 
solution schemes are developed for these problems. The chapter is divided 
into four main sections: (1) general governing equations, (2) governing 
equations for transpiration cooling, (3) governing equations for ablative 
composite plates, and (4) governing equations for restrained thermal growth 
(RTG) testing. For each case, the general governing equations are simplified 
appropriately. The implementation of the solution scheme for each of the 
three problems is also discussed. 

4.1 General Governing Equations 

The general governing equations used here are based on the model of 
McManus [12]. The general governing equations are developed based on two 
models: (1) thermochemical and (2) mechanical models. In the 
thermochemical model, McManus used a control volume approach to obtain 
the local mass and energy balance equations. These two equations were used 
to obtain the internal pressure and temperature distributions inside the 
control volume. The mechanical model was based on the poroelasticity theory 
developed by Biot et. al. [19]. This theory takes into account the internal 
pressure in the calculation of stresses. 
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The unit control volume of a porous solid used by McManus [12] is 
shown in Figure 4-1. In general, a porous solid contains both porous virgin 
and porous charred solids. The volume occupied by the virgin and charred 
solids are denoted by v s and v c , respectively. Initially, the unit control 
volume consists of porous virgin solid and absorbed moisture only. Due to 
heating, some of the porous virgin solid is converted to porous charred solid. 
Gases are then released from the virgin solid into the pores of the virgin and 
charred solid. Moreover, the gases flow through the walls of the unit control 
volume. By applying the principles of conservation of mass and energy to the 
unit control volume, the general 3-D governing equations for pressure and 
temperature are obtained [12]. Then poroelasticity theory is applied to the 
unit control volume to obtain the stresses in the porous solid. In this work, 
only the 1-D form of the general pressure, temperature, and stress governing 
equations are developed. The reason for this simplification is discussed in the 
following section. 

4, LI Ge ome tric -Considerations 

In the first two cases (transpiration cooling and ablative composite 
plate), the structure considered is a thin plate depicted in Figure 4-2. The in- 
plane dimensions are assumed to be much greater than the through- 
thickness dimension. In addition, the boundary conditions, such as 
temperature, pressure, and surface tractions, are uniform over the in-plane 
dimensions. It can then be assumed that all of the derivatives with respect to 
the in-plane dimensions are zero, hence the temperature, pressure, and stress 
vary only in the thickness direction. Therefore, only the 1-D governing 
equations for temperature, pressure, and stress need to be developed. The 
RTG problem can be simplified to 1-D as well. This is 
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Control Volume Control Volume 


Figure 4- 1 . Illustration of the unit control volume. 
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Figure 4-2. Thin plate geometry used in transpiration cooling and ablative 
composite plate where the length (L) and the width (HO are 
much greater than the thickness ( h ). 
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achieved by considering the variations in temperature, pressure, and stresses 
along a strip of the material. More details will be discussed in the later 
sections where the governing equations for the RTG problem are developed. 


The 1-D gas mass balance (continuity) equation is: 

dm g dm g 

ir = - dz +r * 


(4-1) 


where m g is the mass of gas per unit control volume, m g is the mass flux of 
gas, and r g is the rate of gas mass generation per unit control volume. The 
first term of Eqn. 4-1 represents gas mass storage, the second term 
represents the change of mass due to gas flow, and the last term represents 
the generation of gas due to chemical reactions. 

The velocity of the gas flowing inside the porous solid is assumed to 
obey Darcy's Law [31]: 


0 M dz 


(4-2) 


where V D is the area-average gas velocity, y is the permeability of the porous 
solid, fi is the average viscosity of the gas [32], and P is the internal 


pressure. Then the mass flux of the gas is: 


m g =p t V D (4-3) 

where p g is the intrinsic density of gas defined as the density of the gas 
within the pores. The mass of gas per unit control volume is related to its 


intrinsic density by: 

m g = tP* 

where <p is the porosity of the porous solid. 


(4-4) 
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The gas inside the porous solid is assumed to behave ideally, so the 
pressure is given by the ideal gas law as: 

RTm. 


P = 


M<p 


(4-5) 


where R is the universal gas constant, T is the absolute temperature, and M 
is the average molecular weight of the gas [32]. 

The 1-D energy balance equation is: 


^ = - ^( 0 -| :( 0 + Z£ tOTt (4-6) 

where E is the internal energy inside the control volume, q cond is the heat flux 
due to conduction, q conv is the heat flux due to convection, and E g is the heat 
generated by the fcth chemical reaction. Equation 4-6 represents the rate of 
internal energy change inside the control volume due to heat flow by 
conduction and convection and heat generation by chemical reactions. More 
specifically, the terms in Eqn. 4-6 are: 


BE 

dt 



Qcond ^ i ^ 


Qconv = 


I^=S«.a 


(4-7) 


where m i is the mass per unit control volume of the ith substance, /», is the 
enthalpy per unit mass of the ith substance, K z is the area-average thermal 
conductivity in the z direction, h g is the enthalpy per unit mass of the gas, R t 
is the rate of reaction of the Jfcth reaction, and Q is the chemical heat of 
reaction of the kth reaction. By substituting Eqn. 4-7 into Eqn. 4-6 and 
applying the ideal gas assumption, the energy equation (Eqn. 4-6) becomes 
[ 12 ]: 
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Xc, 


P m , 


£L = A(r dT 

dt dz\ 


. BT 


K '-Tz\~ Cp > 1 ”* 17 + S ^2* - 1 r A 


<9z 


(4-8) 


where r is the rate of mass generation of the ith substance, C f is the specific 
heat of the ith substance, and C P is the specific heat of the gas. 

Note that in Eqn. 4-8, the rate of reaction of the ith reaction, R k , is in 
general a function of temperature, internal pressure, state of chemical 
reactions, etc. A reaction rate law is used to predict R t . The rate of mass 
generation of the ith substance in Eqn. 4-8, r t , is given by the sum of 
contributions from all reactions that produced the ith substance 
( r i =1/,^ ). Finally, the char volume, v c , is defined in the unit control 
volume as the ratio of the current mass of charred solid to the mass when all 
reactions are complete. 

Equations 4-1 through 4-8, along with appropriate initial and 
boundary conditions, and reaction laws (discussed later) provide the 
necessary relations to obtain the temperature and internal pressure 
distributions. 


4.1.3 Stress Equations 

In this section, the summation convention over repeated indices is 
assumed and comma indicates spatial derivative. The equation of motion for 
the porous solid without body force is : 


B 2 u t du 


(4-9) 


where p s is the density of the porous solid, u, is the displacement vector of the 
porous solid in the ith direction, c is the damping coefficient, and ( 7$ is the 
total stress tensor. The tr ans ient terms on the left hand side of Eqn. 4-9 are 
due to the inertial and d am p in g effects of the porous solid. The term on the 
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right hand side of Eqn. 4-9 is associated with the forces developed due to the 
deformations of the porous solid. 

The damping term is introduced in order for the solution of Eqn. 4-9 to 
reach steady-state. By nondimensionalizing Eqn. 4-9 (see Appendix A.l), a 
nondimensional damping coefficient is derived [33]: 



where g is the nondimensional damping coefficient, and co n is the circular 
natural frequency and is given by ^E/(p s h 2 ). In the transpiration cooling 
calculation, it will be assumed that the solid response will be underdamped 
( g<l ) and a g value of 0.1 is used. The value of 0.1 for g is unrealisticly 
large. However, in the problems considered the transient effects in stress are 
not very important, hence a large value of g is used to damp out the transient 
effects quickly. 

The total stress tensor, a tj , is defined in reference [12] as: 

(4-11) 

where &* is the mechanical stress tensor and 8 {j is the Kronecker delta. The 
total stress tensor includes contributions from both the mechanical stresses 
and internal pressures. 

The strain-displacement relation is: 

e *=|( M u + “;..) (4 ' 12) 

Strain may be introduced into the porous solid by mechanical stresses, 
internal pressures, temperature, moisture, and chemical (charring) reactions. 
Accordingly, the total strain is: 

e,y = S iJu aZ + A, y A P + a^T + J8 # A (MC) + **Av c (4-13) 
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In Eqn. 4-13, a, /3, and % are th e temperature, moisture, and charring 
expansion coefficients, respectively, and AP, AT, A (MC), and Av c are the 
changes from a reference value of the pressure, temperature, absorbed 
moisture content, and char volume, respectively. The tensor S ijkl is the 
compliance of the porous solid under mechanical loading and A, y is the 
compliance of the porous solid when subjected to an internal pressure. The 
values of these compliance tensors must be determined from experiments. 

Once the temperature and internal pressure distributions (Eqns. 4-1 
through 4-8) are determined, Eqns. 4-9 through 4-13, together with 
appropriate initial and boundary conditions, provide a set of relations 
necessary to obtain the stress distributions inside the porous solid. The 
specific boundary conditions used for each of the three cases will be discussed 
in the following sections. 

4.2 Transpiration Cooling 

To demonstrate the feasibility of using the explicit finite difference 
method (EFDM) to solve a system of partial differential equations on a CM-5, 
a sample problem is first solved. The sample problem chosen is the 
transpiration cooling problem. This problem resembles the ablative 
composite plate problem with the absence of chemical reactions. 

In the transpiration cooling problem, a porous solid, here considered as 
a plate with a through-thickness temperature gradient, is cooled by sending a 
gas flow from the cool side to the hot side. This sort of cooling has found 
applications in the cooling of turbine blades and is under consideration for 
the surfaces of hypersonic vehicles. The problem is illustrated in Figure 4-3. 
The pressures and temperatures are fixed at the boundaries. On one end of 
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the plate U = h), the value of the mechanical stress is fixed (set to 0). On the 
other end (z = 0), the displacement is fixed (set to 0). It is desired in this 
problem to find the pressure, temperature, and stress distributions through 
the thickness of the porous solid matrix. In this problem, all properties 
associated with both the porous solid and cooling gas are assumed to be 
constant. Moreover, the porous solid is isotropic and consists of the porous 
virgin solid only. 


4.2.1 Transpiration Cooling Governing Equations 
In order to obtain the through-thickness temperature, pressure, and 
stress distributions, the general 1-D governing equations need to be 
simplified appropriately for the transpiration cooling problem. Since there 
are no chemical reactions in the transpiration cooling problem, the gas 
generation term in the 1-D continuity equation (Eqn. 4-1) is set to zero: 



(4-14) 


It is assumed that the system consists of two species: an isotropic porous solid 
and a cooling gas. The two energy generation terms due to chemical 
reactions in the energy balance equation (Eqn. 4-8) are set to zero. The 
energy equation (Eqn. 4-8) then simplifies to: 



(4-15) 


where C p is the specific heat of the porous solid and m s is the mass of the 
porous solid per unit control volume. After some algebraic manipulations 
(using Eqns. 4-2 through 4-5) Eqns. 4-14 and 4-15 become: 


i£i = J!L 

dt tpflM 
9T 


^T + 3 


dz 


•r. 


dz 2 


3p >' * A * A''* - * ! 
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, (4-16) 


* (c^l-fJp. + C^pJjL *<** ^ 

The derivation of Eqns. 4-16 is outlined in more details in Appendix. A.2. 
Equations 4-16 are used in the actual computation. 

Simplifying the equations of motion to 1-D: 


dPu, du , 

dt 

d 2 u, du „ 

dt 

d 2 u, du. 


K/ W r VM. 




(4-17) 


P -* f + C A =<T - 


Note that it is assumed in Eqns. 4-17 implicitly that the damping coefficient 
is the same in all three directions. The equations of motion can be simplified 
further by making use of Eqns. 4-11 through 4-13. To be consistent with the 
assumption that the system consists of only the isotropic porous virgin solid 
and the cooling gas, A(MC) and Av e are also set equal to zero in Eqn. 4-13. 
After some algebraic manipulations, Eqns. 4-17 become: 


d 2 u du _ 
n ■■ ■ + c - ■ - * = Du 

Pt ^ dt x 
d 2 u du 

p ’-d?+ c -£~ Du 

P 'l? +c lt° 4 - {A+B)AP ‘- CAT ‘) 


y-a 


(4-18) 


where 
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(4-19) 


^ _ 2u 2 + u- 1 

f(w-l) 

B _ A(u-M) 

{v-l) 
g(v + l) 

(«-l) 

d =_L_ 

2(1 + t>) 

E is the Young's modulus of the porous solid, v is the Poisson's ratio, and the 
pressure compliance tensor, A , is taken from reference [34] for the isotropic 
and dilute porosity case as: 


A = - 


l-2v 

E 


(4-20) 


The derivation of Eqns. 4-18 are outlined in more detailed in Appendix 
A.3. Equations 4-16 through 4-20 are incorporated into a computer code on 
the CM-5 to obtain the through-thickness temperature, pressure, and stress 
distributions. The computer implementation is discussed in the next section. 


4.2.2 CM-5 Implementation - Transpiration Cooling 
In order to obtain the through-thickness temperature and pressure 
distributions, Equations 4-16 are cast into an explicit finite difference form. 
Forward difference is used for the temporal derivation, backward difference, 
with respect to the gas flow direction, is applied to all first order spatial 
derivatives, and central difference is applied to all second order spatial 
derivatives. The finite difference forms are: 
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where p g ' is the intrinsic gas density at the ith node and the ;th time step, 
Tf is the temperature at the /the node and the jth time step, Az is the 
distance between two neighboring nodes, and A t is the time step. These 
approximation schemes are chosen to ensure stability of the explicit finite 
difference method [35]. The values of temperature (T) and intrinsic gas 
density (p g ) are obtained by marching forward in time. The values of the 
internal pressure are calculated using the ideal gas law (Eqn. 4-5) with the 
known value of T and p g at each time step. 

To obtain the through-thickness stress distributions, Eqns. 4-18 are 
cast into an explicit finite difference form. Central difference is applied to 
both the temporal and spatial derivatives for stability reasons [35]: 
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where u c ' is the value of the displacement of the porous solid in the x 
direction at the ith node and the ;th time step, and similarly for u v ' and uj. 
Once the displacements are obtained, the strain-displacement relation (Eqn. 
4-12) is used to obtain the strains. By substituting the strain results into the 
constitutive relations (Eqn. 4-13 with A (MC) and Av f set to zero), the 
through-thickness mechanical stress distributions, <7 % , are obtained. 

The spatial derivatives are computed simultaneously at all nodes. As 
discussed before in Chapter 2, this is achieved by declaring multiple 
conformable arrays for T, p f , and u's in combination with the EOSHIFT 
command [1] to get arrays that contain the appropriate data elements. By 
adding and/or subtracting the arrays according to Eqns. 4-21 and 4-22, the 
finite difference approximation for the spatial derivative are obtained in 
unison for all nodes. This is repeated for each time step until the values of T, 
p s , and u's reach steady state. The results are shown and discussed in 
Chapter 5. 

4.3 Ablative Composite Plate 

To demonstrate that complex physics can be incorporated into the 
EFDM easily, the ablative composite plate problem is solved. Two major 
complexities arise in the ablative composite plate problem that are not 
encountered in the transpiration cooling problem. The first complexity is 
that the material properties are no longer constants. In general, they are 
functions of temperature (D and char volume (v c ) 1121. The second 
complexity is that chemical reactions take place in the ablative composite 
problem. The chemical reaction model used in this study is the two-step 
reaction model by McManus [12]. These additional complexities are 
incorporated into the current model. 


The geometry of a typical ablative composite plate is shown in Figure 
4-4. Composite plies with fibers at an angle 0 (usually 45*) are built up at an 
angle <P to the heated surface. When the surface of the plate is exposed to 
high temperatures (e.g. during rocket firing), the heat conducts into the 
material, causing it to degrade and release gases. These gases create internal 
pressure which can eventually exceed the cross-ply strength of the composite, 
resulting in delaminations (ply-lifts) [12]. 

The boundary conditions are described by referring to Figure 4-4. ' On 
one side of the plate (z = 0) the displacements are fixed. That side is also 
insulated and impermeable to heat and mass flux. The other side ( z = h ) is 
open to ambient environment where the heat flux, pressure, and applied 
mechanical loads are specified as functions of time. Initial conditions are 
specified uniformly through the thickness. Geometry, initial conditions, and 
all of the surface conditions (heat flux, pressure, and applied mechanical 
load) are uniform in the x and y directions. Hence, the solutions (pressure, 
temperature, stress) vary only in the z direction. 


4-3.1 Governing Equations - Ablative Composite. Elate 
There two sources from which gases can be generated in this problem: 
(1) evaporation of absorbed moisture and (2) charring of the porous solid. The 
gas generation due to the first and second source will be denoted by r, and r p , 
respectively. The 1-D continuity equation (Eqn. 4-1) now becomes: 




+ r, + r. 


(4-23) 


Note that the gas generation term in Eqn. 4-1, r t , is now given by r t and r p , 
which denote rates of gas generation due to evaporation and pyrolysis 



Figure 4-4. 


Ablative composite plate geometry. 


reactions, respectively. Moreover, the terms r e and r p are given by the 
following expressions [12]: 


r t = R»m,+(MC)p,R c 

r p=K(p,-p c ) 


(4*24) 


where is the vapor formation rate, R c is the char formation rate, p c is the 
intrinsic density of the charred solid, and A/C is the instantaneous moisture 
content defined as the ratio of the liquid mass over the solid mass. Assuming 
evaporation and charring reactions are temperature-rate dependent, R w and 


R c are given as [12]: 
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where 

is the temperature at which the evaporation begins, T„ 

is the 


temperature at which the evaporation ends, is the temperature at which 
charring begins, and T n is the temperature at which charring ends. Note 
that T^, T„, and T n are functions of pressure. By applying Darcy's law 
(Eqn. 4-2) and the ideal gas law (Eqn. 4-5), the 1-D continuity equation (Eqn. 
4-23) becomes: 
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Eqn. 4-27 is used to obtain the internal pressure distribution. The detailed 
derivation of Eqn. 4-27 is outlined in Appendix B.l. 

The ablative composite plate is a system containing four components: 
(1) porous virgin solid, (2) porous charred solid, (3) absorbed moisture, and (4) 
flowing gases due to evaporation and charring. The sum of the internal 
energies of these four components contribute to the total internal energy of 
the system. Specifying the energy equation (Eqn. 4-8) to this case: 
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where Q. and Q. are the effective heat of charring (pyrolysis) reaction and 
evaporation reaction, respectively. They are given as [12, 22]: 


a - (a +pa -pa-(mc)pa*(mc)pa-(p, -p,)k) 

Q w =(Q w -m J h t + m J h 1 ) 

where Q. is the heat of charring reaction, Q» is the heat of evaporation 
reaction, h } is the specific enthalpy of the virgin solid, h c is the specific 
enthalpy of the charred solid, h p is the specific enthalpy of the pyrolysis gas, 
h, is the specific enthalpy of the absorbed moisture, and h, is the specific 
enthalpy of the evaporation gas. The combined heat capacity of the system, 
7 ] , is given as: 


tj * C,m, + C p m e + C Fi m, + C f m t (4-30) 

where C P is the specific heat of the virgin solid, C p is the specific heat of the 
charred solid, C p is the specific heat of the absorbed moisture, C P/ is the 
specific heat of the flowing gases, m, is the porous virgin solid mass per unit 
control volume, m c is the porous charred solid mass per unit control volume, 
m, is the absorbed moisture mass per unit control volume, and m f is the 
flowing gases mass per unit control volume. Eqn. 4-28 is used to obtain the 



through-thickness temperature distributions of the plate. The derivation of 
Eqn. 4-28 is described in detail in Appendix B.2. 

To obtain the through-thickness stress distributions, it will be shown 
that the time scale associated with the solid deformations is much smaller 
than the time scales associated with the temperature and pressure responses. 
Therefore, the response of the solid matrix is approximately steady-state on 
the time scales of the temperature and pressure responses. This is justified 
by the results of the transpiration cooling problem and will be discussed 
further in Chapter 5. 

The development here for the stress governing equations is based on 
the work by McManus [121. The equations of motion for the solid matrix in 
steady-state condition (Eqn. 4-9) become: 

<7^=0 (4-31) 

which are just the equilibrium equations. As discussed before in section 4.1, 
the problem allows the simplification to 1-D. With this simplification, Eqn. 4- 
31 becomes: 


<y* { =0 (4-32) 

By substituting the definition of the total stress tensor (Eqn. 4-11) into Eqn. 
4-32, the following relations are obtained: 

&* =0 

<7^=0 (4-33) 

For the ablative composite plate problem, the boundary condition at the fixed 
surface ( z = 0 in Figure 4-3) is: 

«,= 0 


(4-34) 


At the free surface ( z = h in Figure 4-3), the boundary condition is: 


P = P„ 


(4-35) 


where f* is the mechanically applied traction on the free surface in the ith 
direction and P b is just the ambient pressure. The results, after integrating 
Eqn. 4-33 through the thickness (in the z direction) and applying the surface 
boundary conditions (Eqn. 4-35), are: 
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The strain-displacement relations for the present problem are: 
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Substituting the strain-displacement relations (Eqn. 4-37) into the 
constitutive equations (Eqn. 4-13) yields: 

0 = + A a AP + a a AT + P a A{MC) + * c Av f 

0 * + A „AP + a„AI + p„A {MC) + Z„*v e 

^ = S^aZ + A a A P + a a AT + P a A{MC) + * a Av f 

^■ = S ydd <rZ+ A n AP + a n AT + P n A(MC) + * n Av c (4-38) 

% = 5 atf o2+A a AP + a a AT+ )3 a A(Af C) + ** Av e 

</Z 

o = + A^AP + c^AT + P^AiMC) + **Av e 

Equations 4-38 are a set of six equations in the six unknowns, the 
three displacements u x , u r , and u z , and the three stresses <T^, <7^, and cx£. 



The other three stresses are given by Eqn. 4-36. Equations 4-38 are solved 
with the boundary conditions specified in Eqns. 4-34. 

4.3.2 CM-5 Implementation - Hybrid Algorithm 

In this section the numerical method used to solve the governing 
equations for pressure and temperature (Eqns. 4-27 and 4-28) is discussed. 
Due to the stringent stability criterion of the explicit finite difference method 
(EFDM), the critical time scale for stability needs to be identified. A method 
for identifying the time scale is proposed. Based on this method, the time 
scales associated with the temperature and pressure governing equations are 
identified. From these time scales, the most critical one for stability is 
determined numerically. A time step is then chosen for the calculation. The 
basic idea of the hybrid algorithm is that implicit finite difference method 
(IFDM) is applied for the regions where the time step is greater than the 
critical time scale, and EFDM is used in the regions where the time step is 
smaller than the critical time scale. How the critical time scale for the 
ablative composite problem is identified is discussed in the following 
paragraphs. 

To identify the time scales, the following two generic equations are 
considered: 

in - v ia 

* ** (4-39) 

du d 2 u 

where u is some scalar quantity, and v and c are constants. Critical time 
scales associated with Eqns. 4-39 are [61: 
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where t cl is the critical time scale for the first equation of Eqns. 4-39, t c2 is the 
critical time scale for the second equation of Eqns. 4-39, and /, is some 
reference length such as the spacing between nodes. 

By comparing the first three terms on the right hand side of the 
pressure governing equation (Eqn. 4-27) to Eqns. 4-39, three time scales' are 
identified by comparison to Eqns. 4-40: 


'pci - 


f Pc2 = 


f = — 1 

VJ - I n. 


<p dz\Tfi 

K 

y_dP_ 
<PH dz 

JL 

Pr I 


(4-41) 


In Eqns. 4-41, t^, t^, and 3 are the time scales associated with the 
governing equation for pressure. The first two equations in Eqns. 4-41 are 
found by comparing the first two terms on the right hand side in the 
governing equation for pressure (Eqn. 4-27) to the first equation in Eqns. 4- 
40. The third equation in Eqns. 4-41 is found by comparing the third term on 
the right hand side of Eqn. 4-27 to the second equation in Eqns. 4-40. The 
time scales associated with the governing equation for temperature (Eqn. 4- 
28) are determined in a similar way: 
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The pressure and temperature governing equations (Eqn. 4-27 and 
Eqn. 4-28) are first solved using the EFDM. The explicit finite difference 
form of the pressure governing equation used in computation is: 
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Equation 4-43 is determined using forward difference for time, backward 
difference, with respect to the gas flow direction, for the spatial gradients of 
the coefficients, forward difference, with respect to the gas flow direction, for 
the pressure gradient, and central difference for the second spatial derivative 
of pressure. This scheme was chosen to maintain stability and accuracy [26, 
36]. The explicit finite difference form of the temperature governing equation 
used in computation is: 
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Equation 4-44 is determined by forward difference for time, backward 
difference, with respect to the gas flow direction, for the spatial gradients of 
the coefficients, forward difference, with respect to the gas flow direction, for 
the spatial gradients of temperature and pressure, and central difference for 
the second spatial derivative of temperature. Again, the scheme is chosen for 
stability reason [26, 36]. 

In the computation of pressure and temperature (Eqns. 4-43 and 4-44), 
the gas mass, m g , is determined from the ideal gas law (Eqn. 4-5). The 
pressures and temperatures are determined by Eqns. 4-43 and 4-44. The six 
time scales associated with the problem are calculated as well. In the 
computation of the time scales, the reference length, /,, is set equal to the 
spacing between the nodes, Az. From the six identified time scales (Eqns. 4- 
41 and. 4-42), it is found in practice that the third equation in Eqn. 4-41 gives 
the smallest time scale [37]. Therefore, pressure response controls the 
critical time scale. This implies that for stability reasons in the EFDM, the 
time step taken in Eqn. 4-43 has to be smaller than the smallest value of the 
critical time scale. 

In Figure 4-5, a semi-log plot of the critical time scale vs. position is 
plotted along with the representative time scales associated with heat 
conduction (/^ = Az 2 / 2KJr f) and porous solid deformation * bzJl^E/p , ). 
This plot is generated by solving Eqns. 4-43 and 4 -44 using from an EFDM 
solution at a sim ulation time of 10 seconds. The solution was for a 3cm thick 
plate, finely meshed with 251 nodes. The problem solved is more completely 
described in section 5.2. Note that the abrupt changes in the critical time 
scale curve are not numerical artifacts. These abrupt changes denote the 
boundaries of different regions of material response. It can be seen from 
Figure 4-5 that the smallest value of the critical time scale occurs on the 
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Figure 4-5. Semi-log plot of the conduction time scale, the porous solid 
time scale, and the critical time scale vs. thickness of the 
ablative composite plate. 



surface and is almost of the same order of magnitude as the time scale 
associated with solid deformation. However, near the insulated end, the 
critical time scale is on the same order of magnitude as the one associated 
with heat conduction. The reason for this large variation is due to the 
through-thickness variation in permeability. Permeability varies by a factor 
on the order of 10 6 from the exposed surface to the insulated end. Depending 
on the degree of permeability, the ablative composite plate can be divided into 
three regions of material response in the thickness direction: (1) char region, 
(2) reaction region, and (3) virgin region. This is illustrated in Figure 4-6. 
This can also be seen from Figure 4-5; the char region extends from about 0 
cm to 0.4 cm, the reaction region extends from about 0.4 cm to about 1.0 cm, 
and the virgin region extends from about 1.0 cm to 3.0 cm. 

If the explicit finite difference equations for pressure and temperature 
(Eqns. 4-43 and 4-44) are applied to all three regions, the time step is 
governed by the extent of the char region. Typically, the char region is not a 
critical region where accurate information on pressure and temperature is 
desired [121 . Therefore, to remove the stringent stability criterion, Eqns. 

4-43 and 4-44 are only applied to the reaction and virgin regions (Figure 4-6). 

In the char region, it is also assumed that the gas storage term can be 
neglected [12]. The pressure distribution is then determined from a finite 
difference formulation of Darcy's law [12]: 
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The ideal gas law (Eqn. 4-5) is used to determine m g , and the result is: 


(4-45) 
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Three regions of an ablative composite plate: (1) char regions, 
(2) reactions region, and (3) virgin region. The reaction region 
is comprised of two sub-regions: (A) evaporation region and (B) 
pyrolysis region. 
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Equation 4-46 results in a quadratic equation for the unknown pressure 
value at the /th node. It can be solved progressively from the node closest to 
the surface to the last node in the char region once the mass flux of gas is 
known. The gas mass flux (m # ) can be calculated by asstiming that all gas 
mass generated in this zone flows to the surface immediately. To calculate 
m g , the gas storage term on the left hand side of the continuity equation (Eqn. 
4-23) is set to zero. The continuity equation (Eqn. 4-23) is then integrated 
numerically over the thickness of the char region. 

The energy equation used in the char region is: 
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Note that the energy equation shown in Eqn. 4-47 is in a different form from 
the one used in the reaction and virgin regions (Eqn. 4-28). The reason is 
that m g is obtained in the char region by integrating the continuity equation 
(Eqn. 4-23). In the other two regions, m t is given by the Darcy’s law (Eqn. 4- 
2). To obtain the explicit finite difference form of Eqn. 4-47 the same scheme 
used to obtain Eqn. 4-44 is also used on Eqn. 4-47: 
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Note that in Eqn. 4-48, m g is already known from the numerical integration 
of the steady-state (gas mass storage term is zero) continuity equation (Eqn. 
4-23). 

There is no time scale associated with the pressure calculation in the 
char region (Eqn. 4-46), so stability is no longer a concern. That is not the 
case for the temperature calculation, since there are still time scales 
associated with the energy equation (Eqn. 4-47). It has been identified 
numerically that the most critical time scale associated with the energy 
equation (Eqn. 4-47) is [371: 



n 

This time scale is associated with the heat convection process. Its magnitude 
is high enough to allow the time steps used in Eqns. 4-43 and 4-44 to be used 
in Eqn. 4-48 as well. A flow-chart of the overall algorithm is shown in Figure 
4-7. 

The stress distributions are obtained by solving the system of 
equations in Eqns. 4-38. A standard LU decomposition scheme is adopted to 
solve the system of equations. Since Eqns. 4-38 are derived based on the 
assumption that the response of the porous solid reaches steady-state 
immediately, there is no time scale associated with the problem, i.e. no 
stability problem. As mentioned before, the steady-state assumption is 
justified by the results from the transpiration problem which are shown in 
Chapter 5. Also, as shown in Figure 4-5, the time scale associated with the 
porous solid response is quite small (about 10~® second). The time scales of 
the other aspects of the ablative composite plate problem are typically on the 



Flow chart of the overall algorithm. 
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Figure 4-7. 







order of 10" 3 second. Therefore, it is quite reasonable to assume the porous 
solid response reaches steady state immediately. 

4.4. Restrained Thermal Growth 

Restrained thermal growth (RTG) experiments were performed by 
Stokes [5] and Hubbert [6] on carbon phenolic specimens preconditioned to 
three different mixture contents: wet (8.0% water by weight), as-received 
(3.6% water by weight) and dry (0.27% water by weight). The specimens 
were cylindrical (0.5 in. diameter and 1.0 in. long). These specimens were 
heated uniformly at a constant rate of temperature change. The stress 
required to hold the specimen at a constant longitudinal strain was recorded 
(Figure 4-8). The specimens were fabricated so that the plane of the carbon 
fabric was perpendicular to the longitudinal axis of the specimen. 

Figure 4-9 is a plot of the restraining stress vs. temperature as 
measured by Stokes [5]. Sullivan et. al. [38] divided the specimen response in 
Figure 4-9 into three temperature regions: the thermoelastic, transition, and 
poroelastic regions. The temperature in the thermoelastic region ranges from 
room temperature (297 K) to approximately 450 K. In this region, the 
measured stress is a result of elastic thermal expansion. At some 
temperature, usually above the cure temperature, secondary hydrogen bonds 
break down and the materials softens. The region above this temperature, 
roughly between 450 K to approximately 600 K for carbon-phenolic materials, 
has been identified as the transition region. Two things happen in this 
region: (1) the measured stress decreases considerably due to material 
softening and (2) the specimen begins to undergo pyrolysis and gases begins 
to accumulate in the specimen. The poroelastic region is at temperatures 
above 600K In the poroelastic region the internal pressure becomes large 
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Figure 4-8. Schematic drawing of an experimental setup for RTG testing. 
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Figure 4-9. 


Experimental results of RTG testing [5]. 



and the specimen response is governed by the internal pressure and the 
material's poroeiastic response. 

Since the specimens were heated uniformly during the RTG tests, the 
specimen response does not change along the specimen length. Therefore, the 
RTG specimen response may be simulated by considering only a cross- 
sectional slice of the cylindrical test specimens. In reference [38], the 
longitudinal stress needed to keep the displacement constant was calculated 
by a finite element method using three-noded, constant-strain, triangular 
finite elements. Since the specimen response in the RTG test is axisymmetric 
about the specimen centerline, only one-quarter of the specimen was 
modeled. The finite element mesh used is shown Figure 4-10. In reference 
[38], the specimen response is simulated using a 2-D code. In this work, due 
to the limitations of the current code, a 1-D approximation of the specimen 
response is used. Referring to Figure 4-10, the 1-D code will be applied along 
the strip at y = 0. The current 1-D code cannot be used for quantitative 
comparison since it is not cast in polar coordinates. However, it can be used 
to perform parametric studies where only qualitative comparison is 
necessary. In this study, the effects of pressure-independent and pressure- 
dependent Arrhenius type chemical reaction rate equations are compared to 
experimental data. 

The pressure distribution along the strip is approximated by using the 
same governing equation for pressure calculation as in the ablative composite 
plate case. The only difference in this case is that the z -direction is now 
replaced by the x -direction. The energy equation is not needed because the 
heating rate, <?T/dr, is specified in the RTG tests. Only stress governing 
equations are needed. 
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In the RTG test, the longitudinal displacement was held fixed. 
Therefore, the specimen can be modeled closely by the plane-strain 
assumption. Under the plane-strain assumption, e^, £ re and are zero. 
With these three strains being zero, the strain-displacement relations (Eqn. 
4-12) become: 
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In reference [38], the material of the RTG specimen is assumed to be 
transversely isotropic with the plane of isotropy ( x-y plane) coincident with 
that of the carbon cloth (Figure 4-10). Directly from the constitutive 
equations for transversely isotropic material under the plane-strain 
assumption [39] , the transverse shear stresses o a and a n are zero. Note 
that there is no y variation along the strip in the x -direction (d/dy = 0). Also, 
referring to Figure 4-10, the displacement u y is zero along the strip. With 
these two simplifications, the strain-displacement relations (Eqn. 4-50) 
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Applying the constitutive equations again [39], <7^ can be shown to be zero. 
The equilibrium equation in this case is: 
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Integrating the equilibrium equation over the length of the strip and applying 
the definition of total stress tensor (Eqn. 4-11) and the traction free boundary 
condition at x=r, an expression for <J a is obtained: 


The following two equations are used to solve for u z and : 

+ A„A/> + ct a AT + (3 a A(MC) + *„Av c 
0 = S„ a oZ + + A „*P + a >? M + P„k(MC) + *„Av, 

along with the boundary condition: 


(4-53) 


(4-54) 


« x =0 r = 0 (4-55) 

Eqn. 4-54 is derived by specializing the constitutive equations in Eqn. 4-13 to 
the RTG problem. With <r£, <r£, and <j£ known, the restraining stress, a", 
can be calculated with the following equation derived from Eqn. 4-13: 

<—T-(S^oZ+ S„a"„ + A„A/> + A(MC) + Xa Av,) (4-56) 


4.4,2 Gqvermng JSqHflfian a for RTG - Arrhenius Type Rate Equations 
In modeling the RTG tests, Sullivan et. al [38] modeled the chemical 
reactions by using a 4-step Arrhenius type rate equation: 




RT 


(4-57) 


where c k is the degree of conversion for the fcth reaction, , n k , and are 
the Arrhenius constants and k * 1, 2, 3, 4. However, McManus [40] argued 
that Eqn. 4-57 is not a proper physical model of the chemical reactions, since 
it can generate gases even if the pressure at the point of generation is higher 
than the saturation pressure of the assumed gaseous substance. For this 
reason, McManus [40] developed a pressure-dependent Arrhenius type rate 


equation such that the gas generation rate above the saturation pressure 
proceeds at an arbitrarily small rate. McManus further assumed that the 
mechanism that limits the reaction rate is due to an increase in the 
activation energy ( £ Jt ) with pressure. Eqn. 4-57 is then solved for E % at a 
known value of pressure: 


E+mRT m {P)ln 


& 


(4-58) 


where T sai (P) is the saturation temperature of the gas at pressure P and <*; 4 is 
the arbitrarily small reaction rate, . In computing the value of , a value 
of 0.01 is selected for [40]. It is assumed that in the first two steps of the 
chemical reactions the gas produced is steam. Therefore, T^P) can be 
determined from the steam tables [41]. The resulting E^ for the first two 

steps of the chemical reaction are tabulated in Table 4.1. The latter two steps 
of the 4-step chemical reaction model are assumed to be pyrolysis reactions 
and are independent of pressure. Therefore, the activation energies are 
constant. The values for the activation energy of the pyrolysis reactions are 
tabulated in Table 4.2. The other Arrhenius constants for pressure- 
dependent and pressure-independent models are tabulated in Tables 4.2 and 
4.3. 


4.4.3 CM-5 Impl ementation - RTG 

Both the pressure-dependent and pressure-independent Arrhenius rate 
equations are incorporated into the pressure computation (Eqn. 4-43). More 
details can be seen in Appendix C. The pressure value is then determined for 


Table 4. 1. Pressure-dependent Activation Energy ( £ ^ ) [40]. 


Pressure 

(atm) 

Reaction 1 
MJ/kg mole 

Reaction 2 
MJ/kg mole 

0 

88.76 

117.24 

1 

88.76 

117.24 

2 

90.64 

117.24 

5 

98.24 

117.24 

10 

104.71 

117.24 

20 

112.11 

117.24 

50 

124.13 

119.23 

100 

135.00 

129.67 

200 

147.94 

142.10 

2000 

147.94 

142.10 j 


Table 4.2. Pressure-dependent Arrhenius Constants [40]. 


Model 

Reaction 



A 



Number 

MJ/k 

:g mole 

1/sec 



Pressure- 1 * 1. 

dependent 2 * 4 

Arrhenius 3 211.4 3. 

4 272.1 5. 

*See Table 4.1. 


Table 4.3. Arrhenius Constants for Pressure-independent Model [40]. 
Model Reaction ^ \ n, 



Number 

MJ/kg mole 

1/sec 

• 

Pressure- 

1 

88.76 

1.20xl0 l ° 

3.5 

independent 

2 

117.24 

4.05xl0 9 

6.5 

Arrhenius 

3 

211.4 

3.86xl0 14 

6.5 


4 

272.1 

5.58xl0 13 

3.3 



each time step. The temperature value at each time step is determined from 
the know heating rate, dT/dl . Then, oZ along the strip y = 0 (Figure 4-9) is 
determined from Eqn. 4-56. The resulting aZ calculated based on pressure- 
dependent and pressure-independent Arrhenius type rate equations are 
presented and discussed in Chapter 5. 


CHAPTER 5 


RESULTS AND DISCUSSIONS 


In this chapter, results from the solutions of the three problems 
considered are presented and discussed. The results of the transpiration 
cooling analysis are used to verify the explicit finite difference method 
(EFDM) and to justify an important simplification used in the solutions of the 
ablative composite plate and RTG problems: the steady-state assumption in 
the stress calculations. Results from the ablative composite plate problem 
are used to demonstrate the capability of the hybrid algorithm to incorporate 
complex physics. An improvement in the model made possible by including 
the gas storage term in the pressure calculation is studied parametrically. 
The RTG results are used to study another, different, complex problem. The 
effects of two different chemical reaction models, a pressure-independent 
Arrhenius type rate equation and a pressure-dependent Arrhenius type rate 
equation are explored with this model. Finally, a performance study of the 
hybrid algorithm is presented and discussed. 

5.1 Transpiration Cooling 

In this problem, a plate made of porous material (aluminum) heated on 
one side is cooled by sending a gas flow (H^O vapor) from the cool side to the 
hot side. The properties of the porous plate and the cooling gas used in the 
computation are listed in Table 5.1. The temperature and pressure values 
are fixed on both sides of the porous plate (Figure 4-2). On one side of the 


Table 5.1 Properties for Porous Plate and Cooling Gas. 


Properties 

Symbol 

Units 

Value 

Specific Heat of Porous 
Solid 

c >, 

J/(kgK) 

167 

Specific Heat of Gas 


J/(kgK) 

2000 

Porosity of Solid 

<P 

dimensionless 

0.05 

Area-average Thermal 
Conductivity 

K. 

W/(m K) 

150 

Permeability of Porous 
Solid 

y 

m 2 

1x10' 17 

Gas Viscosity 

p 

N s/m 2 

lxlO- 5 

Molecular Weight of 
Gas 

M 

kg/kmole 

18 

Young's Modulus of 
Porous Solid 

E 

Pa 

70x10* 

Poisson's Ratio of 
Porous Solid 

V 

dimensionless 

0.3 

Damping Coefficient of 
Porous Solid 

c a 

N s/m 4 

1.0x10® 

Intrinsic Density of 
Porous Solid 

P, 

kg/m 3 

2700 

Plate Thickness 

h 

mm 

30 

Thermal Expansion 

a 

UK 

1.2X10-* 


Coefficient of Porous 
Solid 


a Value of the damping coefficient is computed from Gqn. 4-10 by setting the 
nondimensional damping coefficient, g, to 0.1. 






plate ( z = Omm), the displacements are fixed. On the other side of the plate 
(z = 30mm), the surface tractions are prescribed. The values used for the 
temperature, pressure, and surface tractions on the boundaries are listed in 
Table 5.2. Initially, the temperature and pressure distributions in the porous 
plate are uniform and the porous plate has zero displacement and velocity. 
The initial temperature, pressure, displacement, and velocity values are 
given in Table 5.3. 

To verify the explicit finite difference method (EFDM), the computed 
steady-state temperature and pressure distributions are compared with 
analytical solutions for two different values of boundary pressure. The 
pressures used (0.2 MPa and 2.0 MPa) resulted in steady state gas mass 
fluxes of m g =0.034 and 4.90 kg/s m J , respectively. The analytical solutions for 
the steady-state temperature and pressure distributions are given by the 
following three equations [42, 43]: 

(M) 

(5-2) 
or 




Table 5.2 Boundary Conditions 


Pressure (Pa) 
Temperature (K) 
Tractions (Pa) 


z=0 mm 


P, =2x10* or 2xl0 6 
7>273 

T b = unknown 


z * 30 mm 


7>373 
= 0 


Displacement (m) 


a,= unknown 


Table 5.3 Initial Condition Everywhere in Plate. 


Initial Temperature (K) 273 

Initial Pressure (Pa) 1x10 s 

Initial Displacement (m) 0 

Initial (velocity) 


0 


\2hm,Ru , 

tr m * Pl (5 ' 6) 

<5 ' 7> 

The notation used is given in Table 5.1; A is the universal gas constant. The 
steady-state temperature distribution (Eqn. 5-1) can be derived from the 
governing equation of temperature for the transpiration cooling problem 
(Eqn. 4-15) by setting the left hand side to zero and assuming that all 
coefficients on the right hand side are constant. The first steady-state 
pressure distribution (Eqn. 5-2) is derived using Darcy's law (Eqn. 4-2) by 
assuming constant coefficients and uniform through-thickness temperature 
distribution. The second steady-state pressure distribution (Eqn. 5-3) is 
derived in the same manner as in Eqn. 5-2 with a linearly varying 
temperature distribution through the thickness. Equation 5-2 which is exact 
only for a uniform temperature, is used for the relatively high mass flux case 
( m g = 4.90 kg/s m 2 ) which is insensitive to the temperature distribution. 
Equation 5-3 is used for the relatively low mass flux case (m t = 0.034 kg/s m 2 ) 
where the through-thickness temperature distribution can be approximated 
as linear. The computed results are shown in Figures 5-1 and 5-2 along with 
the analytical solutions for temperature and pressure. As shown in Figures 
5-1 and 5-2, the computed results agree fully with the analytical solutions. 
These results lend confidence in the EFDM solution scheme. 

Transient temperature distributions are shown in Figure 5-3. The 
temperature distribution rises from the initial condition to the steady-state 
distribution in approximately one second. The relatively fast response is due 
to the high thermal conductivity material assumed in the computation. The 
time for the temperature distribution to reach steady-state has been checked 


n r* 



• Numerical (mass flux = 0.034 kg/s* m 2 ) 


Analytical (mass flux = 0.034 kg/s*m 2 ) 



z (mm) 


Figure 5-1. Numerical and analytical steady-state temperature 
distributions for two different gas mass fluxes. 
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Figure 5-2. Numerical and analytical steady-state pressure distributions 
for two different gas mass fluxes. 
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Figure 5-3. Transient temperature distributions for four different 
simulation times: (1) 0.001 second, (2) 0.01 second, (3) 0.10 
second, and (4) 1.00 second. 


Q9 


against an available analytical solution [44]. This time is found to be in good 
agreement with the analytical solution. The analytical solution in reference 
[44] is obtained by ignoring the gas mass flux term in Eqn. 4-15. The 
comparison made here is appropriate, but reasonable, since the gas mass flux 
is relatively small in this case ( m / = 0.034 kg/s m 2 ). 

Transient pressure distributions are shown in Figure 5-4. There are 
two things worth noting: (1) there is a small spatial oscillation in the pressure 
distribution early in the analysis (0.001 second) and (2) the pressure 
distribution reaches steady state about 100 times faster than the 
temperature distribution. The small spatial oscillation in pressure is due to 
a sharp temperature rise near one boundary ( z = 30mm) which increases the 
gas pressure before the flow of gas from the other boundary (z » 0 mm) 
reaches there. This sharp rise in temperature can be seen from Figure 5-3 at 
0.001 second. 

Transient stress distributions ( <f ^ ) are shown in Figure 5-5 for four 
different simulation times: (1) 0.001 second, (2) 0.01 second, (3) 0.10 second, 
and (4) 1.00 second. The steady-state gas mass flux value is 0.034kg/s m 2 
which corresponds to P x = 0.2 MPa in Table 5.2. The stress response is 
relatively fast. This conclusion is reached by noting that the stress 
distribution has the same shape as the pressure distribution by about 0.001 
second. This result is predicted for steady-state stresses by poroelasticity 
theory [16-21]. At a very early time (t = 0.0001 sec), a transient oscillation in 
stress can be seen from Figure 5-6. This is caused by an impulsive 
application of the surface pressure (P,). In the actual situation, the 
application of surface pressure is relatively gradual. When the surface 
pressure is applied gradually, the effects of the oscillatory transient stress is 
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Figure 5-4. Transient pressure distributions for four different simulation 
times: (1) 0.001 second, (2) 0.01 second, (3) 0.010 second, and 
(4) 1.00 second. 
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Figure 5-5. 



Transient stress ( cr") distribution for four different simulation 
times: (1) 0.001 second, (2) 0.01 second, (3) 0.10 second, and (4) 
1.00 second. 
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Figure 5-6. Transient stress ( ) dist 
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not significant. In order to mimic the actual situation more closely in the 
analysis, an unrealistically high damping value (g = 0.1) was used to reduce 
the effects of the transient stress. 

The temperature, pressure, and stresses all reach steady state at 
different times. This observation suggests that each physical process 
(temperature, pressure, and stress) has its own characteristic time scale 
(t !rmperarure « In EFDM, if a single time step At is used for all 
calculations, it must be smaller than the characteristic time scale of the 
fastest physical process (in this case, that of stress) to assure overall stability 
of the numerical scheme. The slower physical processes (in this case, 
temperature and pressure) are calculated with time steps smaller than 
required. Therefore, the computations are done many more times than 
necessary. 

A first step made towards alleviating the stability problem in the 
EFDM is to assume that the stress response is quasistatic. The assumption 
is justified since the oscillatory transient stress is neither important nor 
realistic. The quasistatic stress response assumption is in agreement with 
the assumptions made without justification by Kuhlmann [11], McManus 
[12], and Sullivan [13]. 

The same idea can be extended to the pressure calculation as well. 
Notice that in this analysis the pressure distribution reaches steady-state 
about 100 times faster than that of temperature due to the relatively high 
value of permeability assumed for the porous plate. This observation 
suggests that when the permeability of the porous plate is high enough, it can 
be assumed that the pressure response is quasistatic. 
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5.2 Ablative Composite Plate 

In this sample problem, an ablative composite plate is used to model 
the lining of a rocket nozzle exit cone. The in-plane dimensions of the lining 
are much greater than the through-thickness dimension. Also, heat flux, 
pressure, and surface tractions are assumed to uniformly applied over the 
surface of the lining. Hence a 1-D analysis will suffice to capture the 
response of the lining. The properties of the ablative composite plate 
(FM5055) and the flowing gases are taken from reference [12]. These 
properties are listed in Appendix E. The properties given are on-axis (r,-r 2 - 
x,) properties (see Figure 5-7). In order to relate the on-axis properties to the 
off-axis ( x-y-z ) properties, two rotations need to be performed. The first 
rotation is about the jc 3 axis and the second rotation is about the y axis. The 
method used to rotate the properties is described in Appendix F. 

The geometry of the plate is shown in Figure 4-4. The plate has a 
thickness h equal to 3cm. The ply and fiber angles are <t> = 15* and 9 = 45*, 
respectively. On the exposed side of the plate ( z = 3cm), heat flux, ambient 
pressure, and surface tractions are specified. On the insulated side (z = 
Ocm), the heat flux, gas mass flux, and displacements are specified. Initially, 
the temperature and pressure are uniformly distributed through the 
thickness of the plate. The initial values of the temperature and pressure are 
25 *C and 1 atm, respectively. The initial moisture content ( MC a ) is 3 
percent. The simulation time is 105 seconds. In the first 100 seconds, the 
boundary conditions on the surface (z s 3cm) are held constant. In the last 
five seconds, the ambient temperature, the ambient pressure, and the 
convective heat transfer coefficient are reduced over a five second period to 



<!> ply angle 
0 fiber angle 


systems [12]. The on-axis 
» off-axis system is denoted 



25’C, 1 atm, and 90 W/m 2 , respectively. These conditions are used to 
simulate the firing (time < 100 sec) and shutoff (time > 100 sec) of the rocket 
motor. The boundary and initial conditions are listed in Tables 5.4 and 5.5. 

The temperature and pressure distributions are obtained using the 
hybrid algorithm described in Section 4.3.2 of Chapter 4. A mesh of 251 
nodes and a fixed time step of 0.001 second are used. Whenever the critical 
time scale (&t cr = hz 1 /\Pyl<t>n\) becomes less than 0.001 second at a given node, 
the solution scheme for pressure and temperature at that node is switched 
from EFDM to IFDM automatically by the code. 

Temperature distributions at 10, 50, 80 and 104 seconds are shown in 
Figure 5-8. As expected, the temperature is highest at the surface of the 
plate for the first 100 seconds since the surface is exposed to the ambient 
environment of the rocket motor. After shutoff, the value of temperature on 
the surface drops below that on the inside. This is shown by the temperature 
distribution at 104 second. These four temperature distributions are almost 
the same as those obtained by McManus [12]. Note that the previous 
temperature distributions were obtained without the inclusion of the mass 
storage term in Eqn. 4-27. Hence, the results suggest that the mass storage 
term does not have much effect on the through thickness temperature 
distributions. 

Four pressure distributions are shown in Figure 5-9 at 10, 50, 80, and 
104 seconds. The first three pressure distributions reach a peak then drop 
off. These results are different from the results shown in reference [12] which 
did not include the mass storage term. 
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Table 5.4. Boundary Conditions. 


Boundary Conditions 


z = Omm 


z ss 30mm 


Surface Emissivity 
none 

Boltzmann's Constant 
none 

Surface Absorptivity 
none 

Effective Black Body Temp, 
none 

Convective Heat Transfer Coef. 
none 


Ambient Temperature, 
N/A 


Insulated 
q m = 0.0 W/m 3 

Impermeable, 
m g = 0.0 kg/s m 3 


Surface Emissivity 3 , 
e * S.lxlO* 1 

Boltzmann’s Constant 3 , 
a = 5.670x10-® W/m 3 K 4 

Surface Absorptivity 3 , 
a - 4.9xlO l 

Effective Black Body Temp. 3 , 

T, = 3023 K 

Convective Heat Transfer Coef. b,c , 
h cl x 450 W/m 2 K 
h cl = 100 W/m 3 K 

Ambient Temperature b,c , 

T mX = 3023 K 
T„ 2 = 298 K 

Heat Flux, 

= q^ a +q«,v b 

Ambient Pressure®, 

P M = 10 MPa 
P i2 = 0.1 MPa 


Fixed Displacements, Surface Tractions, 

u, = 0 mm T* = 0 MPa 


a The radiative heat flux, q^, is given by the expression: q^ = - a Tf). 

b The convective heat flux, q^, is given by the expression: q m = hJT^ - T m ). 

c For the first 100 seconds the convective heat transfer coefficient, ambient 
temperature, and ambient pressure are given by the value with the subscript 
one. For times greater than 100 seconds, these values are reduced linearly 
over a five second period to the value with the subscript two. 


Table 5.5. Initial Conditions Everywhere in Plate. 


Initial Pressure 

Initial Temperature 

0.1 MPa 

298 K 
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Temperature (°C) 





However, the locations where the maximum pressures occur are 
approximately the same in the results of this study and that of reference [12]. 
Also notice that when the peak pressure value is reached, it remains 
relatively constant in magnitude and simply moves inward as time 
progresses. This observation implies that a continuous chemical reaction 
front is established [40]. In order to capture this phenomena, a relatively fine 
mesh needs to be used. In this case, a mesh containing 251 nodes is used. 
According to reference [40], a mesh of at least 151 nodes is necessary to 
capture this phenomena. At 104 seconds, the ambient pressure is lowered 
from lOMPa to 2MPa. The peak pressure value at 104 seconds is lower than 
that of previous times. However, the difference between the boundary and 
the maximum internal pressure is the largest at this time. It was found in 
reference [12] that this difference is critical to the mechanism that causes 
delaminations (ply-lifts). The value of pressure differential at which ply-lifts 
occur for the current geometry (<D = 15' and 0 = 45’) has been shown in [12] 
to be 6MPa. Then, according to Figure 5-9, ply-lifts are going to occur at 104 
seconds since the pressure differential at that time is 8MPa. This is in 
agreement with reference [12]. 

The maximum pressure differential predicted by the current pressure 
model is compared with the values predicted by the previous model [12]. The 
results are shown in Figure 5-10. Overall, the current pressure model 
predicts a lower maximum pressure value than the previous model. In 
reference [12], it has been found that for smaller ply angle, <I> = 5\ ply-lifts 
can occur before the rocket motor shuts off. This implies that under the same 
conditions as analyzed in the previous study, the prediction by the current 
model may differ from that of the previous model. Also notice that the 
maximum pressure differential value predicted by the current model is 
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Figure 5-10. Maximum pressure difference comparison between current and 
previous models [37]. 





highest in the first 5 seconds, and then settles to a steady lower value. The 
previous model [12] does not predict the same phenomena. It simply rises to 
a peak value and remains at that value until shutoff. 

Moisture loss and char volume as functions of position through the 
thickness of the plate are plotted in Figures 5-11 and 5-12 at 10, 50, 80, and 
104 seconds. Both are expressed non-dimensionally, with 1.0 indicating 
complete charring or moisture loss. Both evaporation and charring reactions 
occur over a very narrow region and proceeds into the thickness of the plate 
with time. Also, the evaporation front precedes the char front. These 
observations are in agreement with the results from McManus's model [12]. 
However, the extent of both the evaporation and char regions are about 0.5 
cm less than that predicted by McManus at the end of the calculations [12]. 

The through-thickness on-axis stress distributions at 104 seconds are 
shown in Figures 5-13 through 5-16. The stress distributions at 104 seconds 
are chosen because the maximum pressure differential occurs approximately 
at this time (see Figure 5-9). For each figure, two stress distributions are 
shown. The solid line distribution is showing only the internal pressure 
contribution to stress. This is achieved by setting the thermal expansion 
coefficient (a) to zero in the computation, hence excluding any thermal stress 
contributions. The dashed line distribution represents the overall 
contributions from both the internal pressure and thermal stresses. 

For the current geometry (<l> = 15*, 0 = 45*), the mechanical shear 
stresses ( <r? and <r" ) are dominated by thermal stresses as shown in 
Figures 5-13 and 5-14. However, their magnitudes are quite small and do not 
exceed the failure strength of the material. Internal pressure and thermal 
stresses contribute approximately equally to the distribution (Figure 5- 
15). Note the magnitude of varies a great deal through the thickness. 
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Figure 5-11. Moisture loss as a function of through-thickness position at 10, 
50, 80, and 104 seconds. 
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Figure 5-15. Mechanical <T? Xi distribution at 104 seconds. 
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This variation resembles the measured restraining stress at different 
temperatures for the restrained thermal growth (RTG) test. These results 
will be compared to the RTG experimental results in the following section. 
The distribution of cr^ is dominated by the internal pressure contribution 
shown in Figure 5-16. This is to be expected for a relatively small value of 
the ply angle <E> = 15° [12]. In this case, the magnitude of is large enough 
to cause delamination (ply-lifts) according to the maximum stress failure 
criterion. 

5.3 Restrained Thermal Growth (RTG) 

In a restrained thermal growth (RTG) test, a cylindrical specimen is 
heated uniformly at a constant rate. The stress required to hold the 
specimen at a constant longitudinal strain is recorded. A schematic RTG test 
setup is shown in Figure 4-7. The measured restraining stresses to keep the 
specimen at a constant longitudinal strain for different temperatures are 
shown in Figure 4-8. In this section, the computed values of cr^ for the 
ablative composite plate problem are compared qualitatively to the 
experimental RTG results. Then the results of a direct numerical RTG 
simulation are shown and discussed. In the simulation two types of chemical 
reaction model are used: (1) pressure-dependent Arrhenius type rate equation 
and (2) pressure-independent Arrhenius type rate equation. 

The results of an RTG experiment are shown side by side with the <7^ J( 
results from the ablative composite plate problem in Figure 5-17 for a 
qualitative comparison. Sullivan [38] has divided the measured material 
response into three regions: (1) thermoelastic, (2) transition, and (3) 
poroelastic. In the thermoelastic region the measured stress is a result of 
elastic thermal expansion. In the transition region the material response 
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changes from thermoelastic to poroelastic. The initial drop in the measured 
stress is due to material softening. As the temperature increases further, 
gases from the pyrolysis reaction begin to accumulate and the effects of the 
internal pressure become apparent. This causes an increase in the 
magnitude of the restraining stress. Finally, in the poroelastic region, the 
internal pressure becomes excessive and the response of the material is 
highly dependent upon the internal pressure, the material's permeability, 
and the material's poroelastic behavior. As one can see qualitatively, the 
computed results compare well with the experimental results. The three 
regions of material response are captured by the hybrid algorithm. The lack 
of quantitative agreement is due to the differences in geometry and boundary 
conditions between the RTG testing and the ablative composite plate. 

In the numerical RTG simulation, the specimens have a radius of 
0.635cm. The material properties used are the same as the ones listed in 
Table 1 of reference [38]. The boundary and initial conditions used in the 
numerical RTG simulation are listed in Tables 5.6 and 5.7, respectively. The 
values of the Arrhenius constants are listed in Appendix C. 

The computed restraining stresses (<y£) based on the pressure- 
dependent and pressure-independent chemical reaction models are shown 
along with the experimental results in Figure 5-18. Note that on the y-axis 
the absolute value of the compressive restraining stress is used. The 
pressure-dependent Arrhenius reaction model captures the trend of the RTG 
experiment better than the pressure-independent Arrhenius reaction model 
which does not capture the trend at all. The quantitative disagreement is 
due to the difference between the model geometry (a strip) and the real RTG 
geometry (a disc). 
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Table 5.6. Boundary Conditions Used for Numerical RTG Simulation 


z=0 cm 

2=0.635 cm 

9 

Impermeable, m g =0 kg/m* s 

Fixed Ambient Pressure, 


P„=0.1 MPa 

Fixed Displacements, u,= 0 cm 

Zero Surface Tractions, 7^=0 MPa 


117 


Table 5.7. Initial Conditions 3 Used in Numerical RTG Simulation 


Initial Pressure (MPa) 

Initial Temperature (K) 

0.1 

298 


a Initial conditions are uniformly distributed. 




Figure 5-18. Restraining stress vs temperature for pressure-independent 
Arrhenius reaction model, pressure-dependent Arrhenius 
reaction model, and measured results [38] . 



The reason why the pressure-dependent Arrhenius model gives better 
results than the pressure independent Arrhenius model can be deduced from 
Figures 5-19 and 5-20, in which the moisture content and internal pressure 
are plotted with respect to temperature. As shown in Figure 5-19, the 
moisture content reaches zero at a lower temperature for the pressure- 
independent Arrhenius model, causing a large amount of vapor generation at 
relatively low temperature. This results in a much higher internal pressure 
(about two orders of magnitude) as shown in Figure 5-20. Therefore, the 
effects of internal pressure on the restraining stress are manifested much 
sooner in the pressure-independent model than in the pressure-dependent 
model. This can be seen in Figure 5-18 where the poroelastic material 
response is seen at a lower temperature. However, the internal pressure does 
not contribute much to the restraining stress since at low temperatures the 
material is relatively rigid. 

5.4 Performance Study 

The results of a performance study of the hybrid algorithm ablative 
code are presented in this section. All the metrics used to measure the 
performance of a parallel code have been discussed in Section 4 of Chapter 2. 
They are used here to determine the performance of the hybrid algorithm on 
the CM-5 massively parallel computer. 

Two different simulation times are used, 0.2 second and 1.0 second. 
The reason for using two simulation times is to compare the purely parallel 
scheme with a hybrid scheme. In the 0.2 second simulation, none of the 
nodes for the mesh sizes used (51, 71, 91, 101, 111, 121, 131, 141, 151, 501, 
and 1001) has gone unstable, so the solution scheme is purely parallel 
(explicit scheme). However, in the 1.0 second simulation, some of the nodes 
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Moisture Content, MC (100%) 



Temperature (K) 


Figure 5-19. Moisture content, MC, vs. temperature for pressure 
independent Arrhenius reaction model and pressure 
dependent Arrhenius reaction model. 
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Figure 5-20. Internal pressure, P, vs temperature for pressure-independent 
Arrhenius reaction model and pressure-independent Arrhenius 
reaction model. 
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in the mesh have become unstable and a serial solution scheme (implicit 
scheme) has to be used. Table 5.8 lists the number of unstable nodes for each 
mesh size used in the 1.0 second simulation. "Unstable nodes" refers to nodes 
for which the solution scheme for pressure has switched from the EFDM to 
the IFDM. The time steps used in both 0.2 and 1.0 second simulations are 
determined by the third equation in Eqn. 4-41: 



As mentioned before in Chapter 4, this time step has been identified to be the 
most critical one in the ablative problem. The values used for pressure, 
permeability, porosity, and viscosity in Eqn. 5-8 are based on the typical 
values that exist in the evaporation zone (Figure 4-6). The reason for using 
these values is to ensure that the parallel solution scheme (explicit scheme), 
which incorporates complex physics, remains stable in the evaporation zone. 
This way, accurate information on pressure can be obtained which is needed 
to predict failure by ply-lifts. These values are listed in Table 5.9. 

The results of the performance study of the 0.2 second simulations are 
listed in Tables 5.10 and 5.11. The results of the performance study for the 
1.0 second simulation are listed in Tables 5.12 and 5.13. Three different 
numbers of CPUs (1, 32, and 64) are used. 1 CPU is used to simulate the 
code as if it was executed on a serial machine. In tables 5.10, the following 
results for the 0.2 second simulation are shown: mesh size, measured clock 
time for a single, 32, and 64 CPUs, speedup for 32 and 64 CPUs, effective 
Amdahl's fractions for 32 and 64 CPUs, and effective parallelization for 32 
and 64 CPUs. In Table 5.11, the following results for the 0.2 second are 
shown: efficiency for 32 and 64 CPUs, and excess time for 32 and 64 CPUs. 


Table 5.8. Number of Unstable Nodes for Different Mesh Sizes After 1.0 
Second of Simulation Time. 


Mesh 

Size 

Number of 

51 

71 

91 

101 

111 

121 

131 

141 

151 

501 

Unstable 

Nodes 

1 

1 

1 

1 

2 

1 

2 

2 

2 

6 
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In Table 5.12, the same results as in Table 5.10 for the 1.0 second simulation 
are shown. In Table 5.13, the same results as in Table 5.11 for the 1.0 second 
simulation are shown. Note that for the 1.0 second simulation run (Tables 
5.12 and 5.13), performance information for the 1001-noded mesh could not 
be obtained. 

The thing to notice is that for a given simulation time, as the mesh size 
increases, the parallel algorithm outperforms the serial algorithm. However, 
when the problem size is small (smaller number of nodes in this case) the 
serial algorithm is actually better. This can be seen from Table 5.10 for 51 
and 71 node mesh sizes where the effective parallelization (p) is actually 
negative. In Table 5.12, for the 51 nodes mesh size, the effective 
parallelization is also negative. Another measure that also points out the 
serial algorithm is better than the parallel one for relatively small problems 
is the speedup, 5. In Table 5.10, for the mesh sizes of 51 and 71 nodes, the 
speedup in actually less than one which implies that the serial algorithm is 
actually faster. The same is also seen from Table 5.12 for the mesh size of 51 
nodes. 

The second thing worth noting is that in this particular case of the 
hybrid algorithm, larger numbers of CPUs do not always lead to better 
performance in terms of speedup, S. The results shown in Tables 5.10 and 
5.12 indicated that using more CPUs actually slows down the computation. 
However, this result may not be true as the problem size or the simulation 
time is increased further. As the problem size or simulation time is 
increased, more of the CPU time is spent on performing parallel computing. 
As more of the CPUs are participating in the actual parallel computing, the 
performance of the hybrid algorithm using 64 CPUs may finally exceed that 
of the 32 CPUs. Hence, the following observation may be made of the results 
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Table 5.9. Typical Values of the Parameters In the Evaporation Zone. 


P (MPa) 

r (m 2 ) 

<t> 

H (kg/m sec) 

10 

5xl0 18 

0.11 

5x10-* 
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Table 5.10. 0.2 Second Simulation Results Summary I. 


Mesh 

Size 

r,., 

(sec) 

^,V= 32 

(sec) 

(sec) 

*32 

*64 



P» 

P64 

51 

1.15 

1.81 

1.80 

0.64 

0.64 

1.59 

1.57 

-0.59 

-0.57 

71 

2.00 

2.44 

2.43 

0.82 

0.82 

1.23 

1.22 

-0.23 

-0.22 

91 

3.20 

3.00 

3.05 

1.07 

1.05 

0.93 

0.95 

0.07 

0.05 

101 

4.00 

3.54 

3.50 

1.13 

1.15 

0.88 

0.86 

0.12 

0.14 

111 

4.97 

4.08 

4.09 

1.22 

1.22 

0.82 

0.82 

0.18 

0.18 

121 

6.14 

4.65 

4.63 

1.32 

1.33 

0.75 

0.75 

0.25 

0.25 

131 

7.53 

5.26 

5.16 

1.43 

1.46 

0.69 

0.68 

0.31 

0.32 

141 

9.17 

5.63 

5.61 

1.63 

1.64 

0.60 

0.61 

0.40 

0.40 

151 

11.1 

6.43 

6.23 

1.73 

1.78 

0.57 

0.55 

0.43 

0.45 

501 

468.6 

53.3 

54.2 

8.79 

8.64 

0.085 

0.10 

0.915 

0.90 

1001 

4113 

205.6 

208.1 

20.0 

19.8 

0.019 

0.036 

0.981 

0.964 
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Table 5.11. 0.2 Second Simulation Results Summary II. 
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Table 5.12. 1.0 Second Simulation Results S umm ary I. 


Mesh 

Size 

T Sml 

(sec) 

T 

1 .V=32 

(sec) 

7”.V=64 

(sec) 

*^32 

*64 

a t32 

®«64 

Pn 

Pm 

51 

3.16 

4.12 

3.97 

0.77 

0.80 

1.31 

1.26 

-0.31 

-0.26 

71 

7.95 

7.18 

6.84 

1.11 

1.16 

0.90 

0.86 

0.10 

0.14 

91 

13.8 

10.8 

11.0 

1.28 

1.25 

0.78 

0.80 

0.22 

0.20 

101 

17.6 

13.0 

12.9 

1.35 

1.36 

0.73 

0.73 

0.27 

0.2.7 

111 

22.4 

15.9 

16.0 

1.40 

1.40 

0.70 

0.71 

0.30 

0.29 

121 

28.3 

18.3 

18.5 

1.55 

1.53 

0.63 

0.65 

0.37 

0.35 

131 

35.6 

24.5 

21.7 

1.45 

1.64 

0.68 

0.60 

0.32 

0.40 

141 

44.5 

24.7 

24.8 

1.80 

1.79 

0.54 

0.55 

0.46 

0.46 

151 

55.3 

28.4 

28.3 

1.95 

1.95 

0.50 

0.50 

0.50 

0.50 

501 

3316 

306 

319 

10.8 

10.4 

0.06 

0.08 

0.94 

0.92 
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Table 5.13. 1.0 Second Simulation Results Summary II. 


Mesh Size 

e }2 

^64 

32 (sec) 

M (sec) 

51 

0.024 

0.012 

4.02 

3.92. 

71 

0.035 

0.018 

6.93 

6.72 

91 

0.040 

0.020 

10.3 

10.8 

101 

0.042 

0.021 

12.5 

12.7 

111 

0.044 

0.022 

15.2 

15.7 

121 

0.048 

0.024 

17.4 

18.0 

131 

0.045 

0.026 

23.4 

21.1 

141 

0.056 

0.028 

23.3 

24.1 

151 

0.061 

0.030 

26.7 

27.5 

501 

0.340 

0.160 

202 

267 
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in Tables 5.10 and 5.12. Given the competing mechanisms (overhead costs 
and actual parallel computation) which determine the net performance, there 
may be an optimum number of CPUs that should be used to give the best 
performance for the hybrid algorithm. 

The third thing to notice is that the effective parallelization, p, and 
efficiency, e , are consistently better for the 1.0 second simulation than the 
0.2 second simulation (see Tables 5.10 through 5.13). This result is to be 
expected. As the simulation time increases, the number of time integration 
operations increases as well. The time integration operations are done in 
parallel. Therefore, it comes as no surprise that the effective parallelization 
and efficiency are consistently better for the 1.0 second simulation than for 
the 0.2 second simulation. Some of the nodes have become unstable in the 1.0 
second simulation and thus the serial scheme (implicit scheme) has to be 
used for some of the calculation. However, the number of unstable nodes is 
much smaller than the total number of nodes in all cases, so most of the 
computation is still done by the parallel scheme. 

For a large problems, the effective parallelization for the code is 
approximately one, but the efficiency is nowhere near that. This result 
indicates that even though the non-paralleiizabie part of the code is a small 
percentage of the whole code, it drastically reduces the efficiency of the 
parallel computation. This result can be attributed to the fact that a large 
number CPUs sit idle while the non-parallelizable part of the code is running. 
This result is worse for the case of 64 CPUs than for the case of 32 CPUs 
since more CPUs sits idle in the 64 CPUs case. Equation 5-9 below gives an 
expression for the efficiency (e) in terms of the numbers of CPU (AO and the 
effective parallelization (p): 


lOI 


1 


(5-9) 


N{l-p) + p 

The above equation is derived by using Eqns. 2-9 through 2-11. Equation 5-9 
is plotted vs p in Figure 5-21 for two different numbers N (32 and 64). As 
one can see from Figure 5-21, for both 32 and 64 CPUs, e drops very quickly 
as p decreases. Also note that e is consistently worse for the N - 64 case 
than the N = 32 case. 
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Efficiency, 



Figure 5-21. 


Plot of efficiency («) vs effectively parallelization (p) for two 
different numbers of CPU (N = 32 and 64). 
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CHAPTER 6 


CONCLUSIONS AND RECOMMENDATIONS 


6.1 Conclusions 

In this work, numerical solution schemes have been developed on a 
massively parallel computer (a Thinking Machines CM-5) to solve three 
problems: (1) transpiration cooling, (2) ablation of a composite plate, and (3) 
restrained thermal growth (RTG). The numerical solution schemes are based 
on two types of finite difference method: (1) the explicit finite difference 
method (EFDM) and (2) the implicit finite difference method (IFDM). 

It has been found that a solution scheme based on EFDM requires 
relatively small time steps for stability but complex physics can be easily 
incorporated into the solution scheme. On the other hand, a solution scheme 
based on IFDM can use relatively large time steps and still maintain stability 
but complex physics are more difficult to incorporate. 

An efficient solution scheme called the hybrid algorithm has been 
developed to combine the strengths and minimize the weaknesses of both 
solution methods as much as possible. The algorithm was developed by 
simplifying the physics of the problems judiciously. This simplification of 
physics was achieved based on observations that different time scales are 
associated with different physical processes. When the time scale of a 
physical process was much faster than the others and when transient effects 
in that process were not important, the response of the process was assumed 
to be quasistatic. This assumption was used to recast the mathematical 
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model of the problem. In all three problems considered, the stress response 
was always quasistatic and its calculation was decoupled from the others; 
pressure response was sometimes quasistatic and sometimes not. The hybrid 
algorithm developed picks regions for EFDM and IFDM in the pressure 
calculations continuously, moving the boundary as the calculation progresses. 
This solution scheme can be conceptually extended to other engineering 
problems where more than one time scale are involved due to different 
physical processes. 

The algorithm was verified by comparing the numerical results to the 
exact solutions for the transpiration cooling problem. From the numerical 
results of the transpiration cooling problem, the time scales associated with 
stress, pressure, and temperature responses were compared and discussed as 
well. The ablation problem demonstrated the importance of advanced 
physics, such as gas storage terms, to solutions. This problem was also used 
to analyze the performance of the algorithm. The results of the RTG problem 
showed that by incorporating enough physics into the algorithm, the complex 
material responses captured during tests can be reproduced numerically as 
well. 

As expected, it was found that the solutions can be obtained in a 
shorter wall-clock time using the CM-5 than a single CPU computer. The 
reduction of wall-clock time is a function of the size and execution time of the 
problem. The reduction of wall-clock time, however, was not linear with 
respect to the number of CPUs. For example, when using 32 CPUs, one 
would expect to obtain the solutions 32 times faster than when using 1 CPU. 
A speedup of only about 20 times was achieved for 32 CPUs. The speedup did 
not increase when 64 CPUs were used to solve the problem. This result 


seemed to indicate that there is an optimum number of CPUs for a given 
problem. 

By using a standard metric, the hybrid algorithm was found to be 
highly (94 percent) parallelized. However, this did not directly translate into 
high efficiency. When the six percent of the algorithm that ran sequentially 
was running, all but one of the CPUs sit idle and this drastically reduced the 
efficiency of the overall performance. 

6.2 Recommendations 

For future work, the following three recommendations are made: 

(a) Extend the analysis to more complex geometries. For example, 
the analysis can be extended to 2-D, 3-D, and axisymmetric 
problems. EFDM makes this extension relatively easy. 

(b) Include more advanced physics such as: an eroding surface, 
stress-dependent permeability, and Sullivan's thermodynamics 
model. Their effects on the calculated response of ablative 
materials can then be studied. 

(c) Parallelize the code more by using more advanced system 
features designed for the CM-5 to obtain higher efficiency. For 
example, in the ablative composite problem, the IFDM 
subroutine used for sequential pressure calculation could take 
advantage of these features to minimize the idle time. 
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APPENDIX A 


EQUATION DERIVATIONS FOR 
TRANSPIRATION COOLING 


A.1 Nondimensional Damping Coefficient (Eqn. 4-10) 

Eqn. 4-10 is an expression for the nondimensional damping coefficient. 
It is derived as follows. We begin with the equations of motion (Eqn. 4-9): 


<? 2 u du., 

p — r- + c — L = a 
Pt dt 2 dt ,J1 


(A-l) 


Based on the thin plate assumption, the only nonzero spatial derivative is 
with respect to the z direction. By using this assumption, Eqn. A-l reduces 
to the following: 


Ps 


d 2 u { du { 
dt 2 +C dt 


(A-2) 


The displacement vector of the porous solid (u,) and the through-thickness 
direction variable (z) are nondimensionalized by the thickness of the plate 
(h). The nondimensionalized displacement vector and the through-thickness 
direction variable are: 


U : 




(A-3) 


The total stress tensor is nondimensionalized by the Young's modulus of the 
porous solid: 



(A-4) 
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where d i: is the nondimensionalized total stress tensor. Substituting Eqns. 
(A-3) and (A-4) into Eqn. (A*2) yields: 


Pjl 2 , ch 2 3u t _ 

E dt 2 E dt~ - 


iZ.i 


(A-5) 


Note that the coefficient, p s h 2 /E , has units of time squared. This coefficient is 
used to nondimensionalized the time variable: 


t =■ 


Jp,h 2 /E 

Substituting Eqn. A-6 into Eqn. A-5 yields: 

du, c du i _ _ 
dP + 4p s Eth l dt ~ a,t ' 1 

Defining the circular natural frequency as [33]: 




Then Eqn. A- 7 can be written as: 


du t c dui _ . 
dP + P,0> H ft" 0 * 1 


(A-6) 


(A-7) 


(A-8) 


(A-9) 


According to reference [33], the nondimensional damping coefficient is 
defined as: 


2 ? = 


P,®« 


(A-10) 


By using the above definition, Eqn. 4-10 is obtained easily. 


A-2 Derivation of Eqn. 4*16 

To derive the first equation in Eqn. 4-16, we make use of the continuity 
equation (Eqn. 4-14): 
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dt dz 

m ? is related to p^ by the following relation: 

= <pp g 

m g is given by Darcy's law (Eqns. 4-2 and 4-3): 

-YP, dP 


(A-12) 


(A- 13) 

jLi vZ 

By substituting Eqns. A-12 and A-13 into Eqn. A-ll, Eqn. A-ll becomes: . 

itei (A.i4) 

dt dzy n az 

P is given by the ideal gas law: 


^pA = ±( YP, dP 
dt <?z^ n dz 


(A-15) 


Eqn. A-16 is obtained by substituting Eqn. A-15 into Eqn. A- 14 for P: 


* dt pMdzK* dz dz 


(A-16) 


By expanding the derivative out on the right hand side of Eqn. A-16, the first 
equation in Eqn. 4-16 is derived: 


3p, yR \jdp,) 2 , . (drY. 

dt ppM A dz J *1 dz A 


(A-ii) 


To derive the second equation in Eqn. 4-16, we make use of the energy 
equation (Eqn. 4-15): 

t ~ ~ \dYT ,, d l T „ , dT 


[ „ _ \dl „ d l _ „ dl 

(C, m, + Cf m,)— = K ,- r - C r m, — 

m s is related to p, by the following relation: 

">,=(!- <P)P, 


(A-18) 


(A-19) 


t or* 


By substituting Eqns. A-12, A*13, and A-19 into Eqn. A-18, the following 
equation is obtained: 


C'.( l-*)p,+C rt *p}?L.K^ + C'' 


YPg dP\dT 
p dz) dz 


(A-20) 


dt '* dz 1 

By making use of the ideal gas law (Eqn. A- 15) for P, the following equation 
is obtained: 


+ -r t § + <*21, 

Then by expanding the derivative out for the second term on the right hand 
side, the second equation in Eqn. 4-16 is obtained: 


dT 


1 






dz dz 


if dr 


+ «l * 




J J 


(A-22) 


A.3 Derivation of Equations of Motion (Eqn. 4-18) 
To derive Eqn. 4-18, we start with Eqn. 4-17: 

d 2 u x du x 

P '^ + C ^ = <7 ~ 


d\ 
dt 2 
d\ 


du 


(A-23) 


du 

Pt dt 2 +C dt 

By using the definition of the total stress tensor (Eqn. 4-11), the three terms 
on the right hand side of Eqn. A-23 can be written in the following forms: 




(A-24) 


i zin 


By substituting Eqn. A-24 into Eqn. A-23, we obtain another form of Eqn. A- 


p £iL +c iiL. ma - 

d 1 u yi du 

p -ir +c i r^- (A - 25 

d 2 u. du. _ _ 

P -1? +C T=^ P ‘ 

For the transpiration cooling problem, the strain-displacement relations are: 


*xz=*vv = £*v =0 

1 du 

£ = t 

a 2 dz 


1 du 

£ = y - 

*2 dz 


(A-26) 


A 


By specializing the constitutive relations (Eqn. 4-13) for the transpiration 
cooling problem, the following constitutive relations are obtained: 

0 = S cr m + S o m +S &" +AAP+aAT 

w an v ix w xryy w yy ^xxyy zz 

0 =5 <r" +5 a” +S <f” +AAP + aAT 

^xxyy^xx rccc v yy ^xxyy zz 


e =S cT + 5 &* + S &” + A AP + aAT 

zz xxyy xx ^xxyy^yy xux v a ■ * ^ 


(A-27) 


0 =2(5 ro -5 cw )cr; iA * 2 

~ 2 (* 5 ot - Sayy) 0 ^ 

Y n ^2(S aa -S xxyy )o% 

where y a and y n are engineering shear strains. For an isotropic solid, S c 
and are given by the following equations: 


5 =— S =— 

E 


(A-28) 


By substituting the strain-displacement relation (Eqn. A-26) into the 
constitutive relations, Eqn. A-27 can be written as: 


I 




0 =S au (?2+S^ y &Z.+S a>y <j2 + AAP + aAT 
0 BS^oZ+S^oZ+S^oZ + AAP+a&T 

^-^oi + S^ + S^oS + AAP+aAr 
0 =2 ( S -- S -K (A-29) 

j*-2(S„S m )o; 

^ = 2(S ra ' 

The first two equations in Eqn. A-29 are used to solve for <r£ and <f ^ . The 
resulting two equations for cr£ and .are then substituted into the third 
equation in Eqn. A-29. The resulting equation is then solved for o£ . The 
last two equations in Eqn. A-29 are used to solve for <t£ and <r£. The three 
equations for cT a , <r£, and are shown below: 

K = Du *z 

< = Du yi (A-30) 

-b*p -csr) 

where the constants A , B, C, and D are defined in Eqn. 4-19. By 
substituting Eqn. A-30 into Eqn. A-25, the desired result (Eqn. 4-18) is 
obtained: 


d 2 u x du x _ 
P,-^f + c-^- = Du, 


^ 9t 

d\ du y 

at 2 +c at 


p,-^?+ c ^r =Du >.u 




(A-31) 
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APPENDIX B 


DERIVATION OF EQUATIONS FOR 
ABLATIVE COMPOSITE PLATE 


B.l Derivation of Governing Equation for Pressure (Eqn. 4-27) 

We begin the derivation of the governing equation for pressure (Eqn. 
4-27) with the 1-D gas continuity equation for the ablative composite plate 
problem (Eqn. 4-23): 

dm dm 

+r ' +r - <B - U 


m is given by the ideal gas law as: 


P<t>M 

ffl 0 = ■ — 

* RT 


(B-2) 


m is given by the Darcy's law as: 


PJ dP 

m = — 

* p dz 


(B-3) 


By substituting Eqns. B-2 and B-3 into Eqn. B-l, the following equation is 
obtained: 


±(££iL\± PjliP +r +r 

dt\ RT ) * u *J ' ' 


Another form of the ideal is gas law is given below: 


P$ RT 


(B-5) 


By substituting Eqn. B-5 into Eqn. B-4 for p g , the following equation is 


obtained: 



d( PQM \_ d (PM y 

<9rl /?T J " dz{ RT n dz ) Tp + ^ 


(B-6) 


M and R are the only constants in Eqn. B-6. By expanding Eqn. B-6 out and 
rearranging the terms, the governing equation for pressure (Eqn. 4-27) is 


¥ = L±(jL}p¥ + JL¥¥+?zi!P+M( r * 

dt <t> dz\.Tfi) dz 0^ dz dz dz 2 <f>M ' 9 *' 


LiL-LilL 

T dt <p dt 


(B-7) 


B.2 Derivation of Governing Equation for Temperature (Eqn. 4-28) 

We begin the derivation of the governing equation for temperature 
(Eqn. 4-28) with the energy equation (Eqn. 4-8): 




(B-8) 


Recall that in the ablative composite plate analysis, the system consists of 
four components: (1) porous virgin solid, (2) porous charred solid, (3) absorbed 
moisture, and (4) flowing gases (evaporation gas and pyrolysis gas). By 
specializing the energy equation to this system, we obtained the following 


equation: 


(c^m, + C p jn c + C p m { + = 

( r A + r c h c + r t h, + r t h, + r f h ? ) 


(B-9) 


r s , r c , r ( , r e , and r are given in reference [12] by the following expressions: 
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r s = ~PsK 
r c = PA 

r P - (Pi ~ Pc)K (b-io) 

r t = R„m s +{MC)p s R c 
r i-~[K m,+(MC)p,R c ] 

Eqns. B-3 , B-5, and B-10 are substituted into Eqn. B-9 to give the governing 
equation for temperature (Eqn. 4-28): 


ST 

dt 




(B.-ll) 


Q. and Q w are given by the following expressions: 


a = (a + PA -pA-(MC)p,h t +(MC)p,h,-{p, - P c )h p ) 
Qw - [Qw ” m A + m A) 


(B-12) 
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APPENDIX C 


IMPLEMENTATION OF FOUR-STEP 
ARRHENIUS REACTION MODEL 


In this appendix, explicit expressions for r t and r p are given. They are 
used in Eqn. 4-43 to compute the pressure distribution in the numerical 
simulation of restrained thermal growth (RTG) testing. 

The mass generation rates of evaporation and pyrolysis gases, r t and 
r p , are given by the following expressions: 


r , = R»m t + (MC)p,R c 

r P = K{p,-Pc) 

R^ and R c are given by the following expressions: 


(C-l) 


(C-2) 


R w = vW + v?R l 
R c = v^R 2 + v^ + viR, 

where (*=1, 2, 3, and 4) are given by Eqn. 4-57 and the values for the 
constants in Eqn. 4-57 are listed in Tables 4.1- 4.3. v k and v c k are given by 
the following expressions: 


m, 




v: = 


(C-3) 


where is the mass of absorbed moisture converted to gas by the 
completion of the fcth reaction, m Ct is the mass of porous charred solid created 
by the completion of the kth reaction, is the total mass of absorbed 
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moisture converted to gas, and m Cj is the total mass of the porous charred 
solid created by all reactions. The values of m, t , m Ci , m, , and m Ca used in the 
RTG numerical simulation are listed in Table C-l [40] . 
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Table C-l Reaction Constants 


Reaction Number 


m c 


kg/m 3 

kg/m 3 

1 

22 

0 

2 

27 

80 

3 

0 

618 

4 

0 

314 


m,' =47 kg/m 3 

m C ' =1012 kg/m 3 
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APPENDIX D 

MATERIAL PROPERTIES OF FM5055 


Table D-l. Material properties of FM5055 carbon-phenolic - permeability 
data for various values of char volume, v e . 


Permeability Units (m 2 x 10' 18 ) 

Char 

Volume, v c 

r* 



0.0 

5 

9 

0.01 

0.1 

6 

13 

0.23 

0.2 

15 

15 

0.90 

0.3 

39 

39 

3.5 

0.4 

112 

112 

13 

0.5 

320 

320 

50 

0.6 

940 

940 

190 

0.7 

2700 

2700 

720 

0.8 

7800 

7800 

2700 

0.9 

23000 

23000 

10000 


1.0 


65000 


65000 


39000 


Table D-2. Material properties of FM5055 carbon phenolic - thermal 
conductivity data for various values of temperature, T. 


Thermal Conductivity Units (W/m 2 °C) 


Temp., 

7TC) 



• 

0 

1.08 

1.08 

0.80 

200 

1.50 

1.50 

0.93 

275 

1.55 

1.55 

1.00 

400 

1.55 

1.55 

1.00 

600 

1.55 

1.55 

1.00 

827 

1.55 

1.55 

1.00 

1227 

2.59 

2.59 

1.60 

1727+ 

3.00 

3.00 

1.60 
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Table D-3. Material properties of FM5055 carbon phenolic - specific heat 

capacity 3 for virgin and charred solid for various values of 
temperature, T . 


Specific Heat Units (J/kg°C) 


Temp., 

TCO 

c >. 


0 

880 

1800 

100 

1165 

1800 

200 

1450 

1800 

300 

1475 

1800 

400 

1500 

1800 

700 

1500 

1800 

800+ 

1500 

1900 


a The specific heat capacities for the flowing gas and absorbed moisture are 
assumed to be constant. In computation, the value used for the flowing gas is 
C P =2000 J/kg°C and the value used for the absorbed moisture is C P{ =4200 

J/kgC. 


Table D-4. Viscosities of flowing gas for two values of temperature. T. 


Viscosities Units (kg/m sec x 10' 8 ) 

Temp., CC) 

Viscosity, (i 

-273 

0.7975 

2727 

8.2975 
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Table D-5. Temperature at which reactions begin and end (°C) and heats of 
reaction (MJ/kg) for various values of pressure. 


Moisture evaporation reaction 

Charring Reaction 

P(MPa) 


T 

* ew 

a 



Qc 

0.1 

100 

150 

*2.251 

400 

538 

-0.234 

0.2 

120 

170 

-2.197 

400 

538 

-0.234 

0.5 

152 

202 

-2.108 

400 

538 

-0.234 

1.0 

180 

230 

-2.008 

400 

538 

-0.234 

2.0 

212 

262 

-1.871 

400 

538 

-0.234 

5.0 

264 

314 

-1.623 

400 

538 

-0.234 

10.0 

311 

361 

-1.351 

400 

538 

-0.234 

20.0+ 

367 

417 

-0.551 

400 

538 

-0.234 
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Table D-6. 

Material Properties of FM5055 - Young's and shear Moduli 
for various values of temperature, T. 

Moduli Units (GPa) 

Temp.CC) 


E , 





23 

17.93 

17.93 

16.55 

6.90 

5.17 

5.17 

93 

17.24 

17.24 

12.41 

6.55 

4.83 

4.83 

204 

13.10 

13.10 

5.52 

3.24 

2.76 

2.76 

315 

8.96 

8.96 

1.38 

2.21 

1.03 

1.03 

426 

6.90 

6.90 

0.55 

1.86 

1.03 

1.03 

537 

6.55 

6.55 

0.34 

1.72 

1.10 

1.10 

815 

7.58 

7.58 

0.34 

1.79 

1.38 

1.38 

1093 

8.27 

8.27 

0.34 

1.93 

1.52 

1.52 

1648 

6.90 

6.90 

0.34 

2.07 

1.52 

1.52 

2204 

2.76 

2.76 

0.34 

1.52 

1.24 

1.24 

2760 

1.38 

1.38 

0.34 

0.69 

0.48 

0.48 
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Table D-7. Material properties of FM5055 - Poisson's ratios for 
various values of temperature, T. 


Temp. (°C) 




23 

0.32 

0.24 

0.24 

93 

0.29 

0.20 

0.20 

204 

0.23 

0.18 

0.18 

315 

0.13 

0.12 

0.12 

426 

0.08 

0.06 

0.06 

537 

0.06 

0.05 

0.05 

815 

0.05 

0.01 

0.01 

1093 

0.06 

0.01 

0.01 

1648 

0.07 

0.01 

0.01 

2204 

0.08 

0.02 

0.02 

2760 

0.09 

0.03 

0.03 
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Table D-8. Material properties of FM5055 - thermal expansions v for 
various values of temperature, T. 


Thermal Expansion Units (m/m x 10’ 3 ) 


TCC ) 

AT('C) 

a AT 

z \ 

a H AT 

a x AT 

23 

0 

0.0 

0.0 

0.0 

93 

60 

0.8 

0.8 

1.0 

204 

181 

1.6 

1.6 

4.0 

315 

292 

2.0 

2.0 

7.0 

426 

403 

2.0 

2.0 

10.0 

537 

514 

1.0 

1.0 

13.0 

815 

792 

0.0 

0.0 

12.0 

1093 

1060 

0.0 

0.0 

-6.0 

1648 

1625 

2.5 

2.5 

-21.0 

2204 

2181 

17.5 

17.5 

-30.0 

2760 

2737 

28.0 

28.0 

-60.0 
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APPENDIX E 

ROTATIONAL MATRICES 


Four matrices are listed in this appendix. The first two (R 2 and R y ) 
are associated with rotating second order tensorial quantities about the z and 
y axes. The latter two (X z and X y ) are associated specifically with rotating 
engineering strains about the z and y axes. Rotation of stress, engineering 
strain, and stiffness and compliance using engineering strain are illustrated 
here. Some relations used in this thesis are expressed in terms of tensor 
strain and tensor stiffness and compliance; the relation between tensor and 
engineering quantities is discussed thoroughly in [39]. 


m 2 n 2 0 0 0 2 mn 

n 2 m 2 0 0 0 -2mn 

0 0 1 0 0 0 

0 0 0 m -n 0 

0 0 0 n m 0 

-mn mn 0 0 0 m 2 -n 2 


(E-l) 


V 2 

0 

*' 2 

0 

-2m n' 

o' 


0 

1 

0 

0 

0 

0 


n' 2 

0 

m' 2 

0 

2mV 

0 

(E-2) 

0 

0 

0 

m' 

0 

n' 


mV 

0 

-m'n 

0 

_'2 >1 
Ttl — 71 

0 


0 

0 

0 

-n' 

0 

m' 
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m 2 n 2 0 0 0 mn 

n 2 m 2 0 0 0 -mn 

0 0 1 ° ° ° (E-3) 

0 0 0 m -n 0 

0 0 0 n m 0 

-2 mn 2 mn 0 0 0 m 2 - n 2 


0 1 0 0 0 0 

0 m- 0 m v 0 

0 0 0 m' 0 n ' 

2 mV 0 -2m V 0 m ,2 -n' 2 0 

0 0 0 -n' 0 m'_ 

where m is short for cosQ, n is short for sin 0 , m' is short for co5<I>, and n is 
short for sin <l> . 

To illustrate how the matrices in Eqns. E-l and E-2 are used to relate 
off-axis second order tensors to on-axis second order tensors, we use the two 
matrices to relate the on-axis (x, -x 2 -x 3 ) stress tensor to the off-axis 
( x - y - z ) stress tensor: 


°^= R A cr »-, 

where o ofr and are given by the following equations: 


= 


a °" “ a. 


Similarly, the on-axis engineering strains can be related to off-axis 
engineering strain by the same expression shown in Eqn. E-5, except now R v 
and R, are replaced by X. and X y , respectively. 
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Eqn. E-7 gives the expression for relating the off-axis compliance 
tensor to the on-axis compliance tensor: 


s* = R,R,s„x;'x; 1 (E-7: 

where XT 1 is the inverse matrix for X 2 and X" 1 is the inverse matrix for X v . 

The inverse matrices for R 2 , R y , X z , and X v are listed below in Eqns. 

E-8 through E-ll: 


'm 2 n 2 0 0 0 -2 mn 

n 2 m 2 0 0 0 2 mn 

0 0 1 0 0 0 

0 0 0 m n 0 

0 0 0 -n m 0 

mn -mn 0 0 0 m 2 -n 2 


' m' 2 0 n' 2 0 

0 10 0 
n’ 2 0 m' 2 0 

0 0 0 m' 

-mri 0 m'n' 0 
0 0 0 n' 


2mV 0 

0 0 

-2 m' ri 0 

0 -n 

m' 2 - n' 2 0 

0 m 


' m 2 n 2 0 0 0 -mn 

n 2 m 2 0 0 0 mn 

0 0 1 0 0 0 

0 0 0 m n 0 

0 0 0 —n m 0 

2 mn —2mn 0 0 0 m 2 —n _ 

m' 2 0 n' 2 0 m'n' 0 

0 1 0 0 0 0 

n' 2 0 m' 2 0 -m'n’ 0 

0 0 0 m' 0 —ri 

-2m' ri 0 2m' ri 0 m' 2 -ri 2 0 

0 0 0 ri 0 m’ 


(E-8) 


(E-9) 


(E-10) 


(E-ll) 
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To relate the off-axis stress tensor to the on-axis stress tensor, the 
following expression is used: 

(E-i2) 

Eqn. E-12 is also valid for all second order tensorial quantities. To relate the 
off-axis engineering strain tensor to the off-axis engineering tensor, replace 
R:‘ and R^ 1 by X7 1 and X^ 1 , respectively The off-axis compliance tensor is 
related to the on-axis tensor by the following equation: 

s 0 „ = r;'r;‘s^x > x j (E- 13 ) 
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APPENDIX F 


INPUT AND OUTPUT FILES OF 
TRANSPIRATION COOLING PROGRAM 


There is one input file and two output files for the transpiration cooling 
program. The name of the input file is "input.txt.” The names of the two 
output files are "outputsl.txt" and "outputs2.txt." The formats and contents 
of these file are shown below. Line numbers are shown for convenience; they 
are not included in the input files. Input values should be valid FORTRAN 
real or integer numbers. Lines are separated by carriage returns, values on 
the same line should be separated by spaces. 

F.l Format of Input File "Input.txt" 

Input.txt 
Line numbers: 

1. CPS 

Specific heat capacity of porous solid (J/kg °C) 

2. CPG 

Specific heat capacity of gas (J/kg 'C) 

3. PORE 

Porosity of porous solid 

4. COND 

Area-average coefficient of thermal conductivity (W/m 2 *C) 


5. PERM 

Permeability of porous solid (m 2 ) 

6. MU 

Gas viscosity (kg/m sec) 

7. MW 

Molecular weight of gas (kg/kmole) 

8. E 

Young's modulus of porous solid (Pa) 

9. v 

Poisson's ratio of porous solid 

10. ALPHA 

Coefficient of thermal expansion of porous solid (lTC) 

11. DAMP 

Damping coefficient (N sec/m 4 ) 

12. MS 

Density of porous solid (kg/m 3 ) 

13. XO 

Thickness of plate (m) 

14. NODE 

Number of nodes in the mesh 

15. - DT 

Time step (sec) 

16. TIME 

Length of simulation time (sec) 

17. TA 

Boundary temperature at z = 0.0 cm (Pa) 
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18. 


TB 

Boundary temperature at z = XO cm (Pa) 

19. PA 

Boundary pressure at z = 0.0 cm (Pa) 

20. PB 

Boundary pressure at z = XO cm (Pa) 

21. TIB 

Surface traction T\ at z = XO cm (Pa) 

22. T2B 

Surface traction T* at z = XO cm (Pa) 

23. T3B 

Surface traction T* at z = XO cm (Pa) 

24. TINIT 

Initial temperature value (*C) 

25. PINIT 

Initial pressure value (Pa) 

F.2 Format of Ouput Files 
Outputsl.txt 


Col. 1 

Col. 2 

Col. 3 

Col. 4 

Col. 5 

Col. 6 

Col. 7 

Position 

Pressure 

Temp. 

Displ. u, 

Displ. Uj 

Displ. Uj 

Sim. 

(m) 

(Pa) 

(K) 

(m) 

(m) 

(m) 

Time (sec) 
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Outputs2.txt 


Col. 1 

Col. 2 

Col. 3 

Col. 4 

Col. 5 

Col. 6 

Col. 7 

Position 




< 


Sim. 

(m) 

(Pa) 

(Pa) 

(Pa) 

(Pa) 

(Pa) 

Time(sec) 
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APPENDIX G 


INPUT AND OUTPUT FILES OF ABLATIVE 
COMPOSITE PLATE PROGRAM 


The input files for the ablative composite plate program consists of six 
categories: (1) general information, (2) mesh information, (3) material 
properties, (4) restart conditions, (5) boundary conditions, and (6) chemical 
reaction constants. In the general information category, there is one input 
file named: "input.in." In the mesh information category, there is one input 
file named: "mesh.in." In the material properties category, there are ten 
input files named: (1) "alpha.in," (2) "cond.in," (3) "denpor.in," (4) "mole.in," 
(5) "poisson.in," (6) "shear.in," (7) "speh.in," (8) "perm.in,"(9) "visc.in," and (10) 
"young.in." In the initial conditions category, there are four input files 
named: (1) "rechar.in," (2) "remois.in," (3) "retpinit.in," and (4) "trate.in." 
These files can also be used to restart runs, and are refered to as restart files. 
In the boundary conditions category, there are three input files named: (1) 
"heatb.flux," (2) "presb.in," and (3) "traction.in." In the chemical reaction 
category, there are two input files named: (1) "ccharl.in," and (2) "wchar2.in." 
There are six output files from the program and they are: (1) "char .out," (2) 
"press. out," (3) "stability.out" (4) "stress. out," (5) "temp. out," and (6) 
"volat.out." The contents and formats of the input and output files are shown 
in the following sections. Line numbers are shown for convienience; they are 
not included in the input files. Input values should be valid FORTRAN real 
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or integer numbers. Lines are separted by carriage returns, values on the 
same line should be separted by spaces. 

G.l General Information Input File 

input.in 
Line numbers: 

1. SEP2 

Number of nodes using IFDM solution scheme for pressure (Eqn. 4-46) 

2. SEPT2 

Number of nodes using 2nd EFDM solution scheme for temperature 
(Eqn. 4-48) 

3. NF 

Number of specified output time 

4. TSPEC(l) 

First specified time (sec) 


TSPEC(NF) 

NFth specified time (sec) 

5. GN 

Number of char volumes at which the permeability of the porous solid 
are specified 

6. KN 

Number of temperatures at which the coefficients of thermal 
conductivity are specified 
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7. NET 

Number of temperatures at which mechanical properties are 
specified 

8. NEC 

Number of char volumes at which mechanical properties are 
specified 

9. CN 

Number of temperatures at which the specific heats of the porous 
virgin and charred solids are specified 

10. MN 

Number of temperatures at which the gas viscosities are specified 

11. TYPEC 

Type of char reaction. For ablative composite plate problem enter 1 

12. TYPEW 

Type of evaporation reaction. For ablative composite plate problem 
enter 2 

13. CHEMNC 

Number of pressures at which the values of Q w > and T„ are 

specified 

14. NODE 

Number of nodes in the mesh 


167 


G.2 Mesh Information Input File 

mesh, in 

1. XO 

Thickness of plate (m) 

2. PT 

Ply angle (degrees) 

3. FT 

Fiber angle (degrees) 

4. DT 

Time step (sec) 

5. TIME 

Length of simulation time (sec) 
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G.3 Material Properties Input Files 

alpha.in (coefficients of thermal expansion of porous solid damaged and 
undamaged (1/°C) ) 

Undamaged Damaged 


Col. 1 

a x,x, ^.V,) 

Col. 2 
a**, (r„v,) 

Col. 3 

«x,x, (^.V,) 

Col. 4 

a x,x, <?>,) 

Col. 5 
Ox^ (T„v,) 

Col. 6 
a XiXy \T liV] ) 

«X,X, 

^NET» V 1^ 

(^NET> V |) 

( ^NET > V 1 ^ 

“X.X, 

(^NET> V |) 

a x,xj 

(^NET> V l) 

«x 5 x, 

(^NET> V l) 

«x,x, 

V NEC^ 

®*»Xj 

(^1» V NEC^ 

(^l» V NEC^ 

(^1> V NEC^ 

«« 

(^I* V NEc) 

«x^ 

Vj 

(^l* V NEC^ 

«x,x, 


«x* 

«x,x, 

«x,x, 

«x* 


(^NET> V NEc) (^NET> V NEC 


» v nec) (Tnct.Vnec) (T'net » V NEC^ ( t net> v nbc'> W mt V„c) 
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cond.in 


(coefficients of thermal conductivity of porous solid 
damaged and undamaged (W/m 2 *C) and temperatures at which 
the coefficients of thermal conductivity are specified CC) ) 
Undamaged Damaged 


Col. 1 

Col. 2 

(7,) 

Col. 3 

(T.) 

Col. 4 

K, (T.) 

Col. 5 

K, (r,) 

Col. 6 
(T.) 

Col. 7 

r, 

(W 

( T^) 


(r„) 

(r„) 


^KN 


170 


denpor.in 


(intrinsic densities of virgin 
of virgin and charred solids, 


and charred solids (kg/m 3 ), porosity 
and initial moisture content) 

Col. 4 


Col. 1 


Col. 2 


Col. 3 


mole. in 


(Molecular weights of evaporation and pyrolysis gases 
(kg/kmole) ) 


Col. 1 
K 


Col. 2 
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poisson. m 


(Poisson’s ratio of porous solid undamaged and damaged) 


Undamaged 


Damaged 
Col. 5 


Col. 6 


Col. 1 Col. 2 Col. 3 Col. 4 

v ( r. , v. ) u^dpv,) U xilj (r„v,) v XtZi (r„v,) (r„v,) v Vj (r„v,) 




u 




u 


Vj 


u 


Vi 


u 


V> 


V 


Vj 


(T-njt.v,) o'nct.v,) (T^.v,) OW.v,) (r^.v,) (7-^.v,) 






(^l> V NEc) ^1 • V 1 


NEC ' 


) (r„V nec) (^.V^) (^*I> Vnec) 


t) 


Vi 


V 


V) 


V) 




Vi 

(W.v^) (r^.v^e) (T^.v^) (?»,»«) ‘Vr.'W) (W.'-nec) 


‘rt 


V 


Vj 


V 


*t *7 


V 


Vj 


Vj 
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shear.in 


(Shear moduli of porous solid undamaged and damaged (Pa) ) 
Undamaged Damaged 

Col. 1 Col. 2 Col. 3 Col. 4 Col. 5 Col. 6 

G Vi (r,,v,) G IzIi ( T t , v, ) G Vi (T x ,v,) G Z]Ij ( 7*1 » v , ) (7>,) (r„v,) 


V! 


V) 


*2*J 


*1*2 


V) 


* 2*3 


(T NCT ,V,) ^.V,) ^NET. V l) (r NET» V |) ^NET.V.) 


* 1*2 

(^..Vnec) (T’lVnec) (T^Vmz) (^.Vnec) 


* 2*3 


* 2*3 


* 1*2 


* 2*3 


* 2*3 


* 1*2 

(Tnet^vnhc) (Tnct.vnec) (r^.v^) (r^.v^) (r^.v^) 


* 2*3 


* 2*3 


* 1*2 


* 2*3 


* 2*3 
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speh.in 


(specific heats for porous virgin and charred solids (J/kg °C), 
temperatures at which the specific heats of porous solid are 
specified CC), specific heats for absorbed moisture, evaporation 
gas, and pyrolysis gas (J/kg °C) ) 


Col. 1 

Col. 2 

Col. 3 

c„, W 

c Pt (V 

7-, 

c„, (Tc*) 

Cp, (r„) 

Tcs 

Q>, 


c,. 
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perm, in 


(permeability of porous solid undamaged and damaged (m 2 ) and 
char volumes at which the permeability of porous solid are 
specified) 

Undamaged Damaged 


Col. 1 

Col. 2 
7 x 2 (v.) 

Col. 3 
7x, <v.) 

Col. 4 
7x, (v,) 

Col. 5 

7xj 

Col. 6 
7x, < v .) 

Col. 7 
v . 

7*. (v GN ) 

7x 5 ( v gn) 

7x, ( V Gn) 

7x, ( v gn) 

7x 2 ( v gn) 

7x, ( v gn) 

V GN 
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(viscosities of evaporation and pyrolysis gas (kg/m sec) and 
temperatures at which the gas viscosities are specified ( C) ) 


Col. 1 

Col. 2 

Col. 3 

M r «) 

m p (t,) 

r, 

Me 

M p 

^MN 
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young. in 


(Young's moduli of porous solid undamaged and damaged (Pa), 
char volumes at which Young's moduli are specified, and 
temperatures at which Young's moduli are specified CC) ) 
Undamaged Damaged 


E x (r„v,) E Xi (r,,v,) E Xi (r„v,) E Xi ( 7 t , v, ) e h ( 7, , v, ) e X} (7,,v,) 


(7"nET> V 1 ) ( ^NCT » V 1 ^ (^NET> V 1^ (^*NET* V 1^ ^NET > V 1 ) ^NET* V l) 


(7, V N£C^ (^1 V NEC^ (^l V,^) (7, V^) (7, V,^) (7, V,^) 


( 7n£j , v nec ^ (^NET> V NEc) ^ ^NET > V NEC ^ (^NET* V NEC^ (7 'nET» V NEc) ( ^NFf ’ V NEC ^ 
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G.4 Restart Conditions Input Files 

rechar.in (restart through-thickness char and solid volume distributions) 

VC(1) 


VC(NODE) 

VS(1) 


VS(NODE) 
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remois.in (restart through-thickness moisture content distribution) 

MC(1) 

MC(NODE) 
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retpinit.in (restart through-thickness temperature CC) and pressure (Pa) 
distributions) 


T(l) 


T(NODE) 

P(l) 


P(NODE) 


trate.in 


(restart through-thickness rate of temperature change idT/dt) 
distribution CC/sec) ) 


TRATE(l) 


TRATE(NODE) 
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G.5 Boundary Conditions Input Files 

heatb.flux (radiative and convective heat transfer parameters at exposed 


surface) 

Col. 1 Col. 2 Col. 3 Col. 4 

a e a T ‘ 

W/m 2 K 4 K 


Col. 5 

Col. 6 

Col. 7 

Col. 8 

Kx 

Ki 

T m 

“l 

T mi 

W/m’K 

W/m^K 

K 

K 
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presb.in 


(pressure at exposed surface and ambient pressure 
(Pa)) 


Col. 1 

P bi 


Col. 2 


b 2 


Col. 3 


amb 
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traction.in (surface tractions at exposed surface (Pa) ) 
TIB 
T2B 
T3B 
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G.6 Chemical Reaction Constants Input Files 

ccharl.in (chemical reaction constants for pyrolysis reaction) 

Col. 1 Col. 2 Col ‘ 3 

r* co r ec co a (J) 
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wchar2.in (chemical reaction constants for evaporation reaction and 
pressures at which these reaction constants are specified) 
Col. 1 Col. 2 Col. 3 Col. 4 

T ^ (P,) CC) T„ (Px) CC) a (Pi) ( J ) Pi (Pa) 

T ^ (PchEMNc) CC) T„ (PcHEMNc) ('C) & (PcHEMNc) (J) PcHEMKC (P») 
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G.7 Output Files 

char.out (through-thickness distributions of virgin and charred solid 
volumes at specified output times (sec) ) 

Col. 1 Col. 2 Col. 3 Col. 4 

Position (m) v, (TSPEC(l) ) v c (TSPEC(l) ) TSPEC(l) 

Position (m) v J (TSPEC(NF) ) v c (TSPEC(NF) ) TSPEC(NF)' 
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press.out (through-thickness distributions of pressure (Pa) at specified 
output times and time step used in simulation (sec) ) 


Col. 1 Col. 2 


Position (m) P (TSPEC(l) ) 


Col. 3 
TSPEC(l) 


Col. 4 
DT 


Position (m) P (TSPEC(NF) ) TSPEC(NF) 


DT 
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stability.out (number of nodes that used Eqn. 4-46 and Eqn. 4-48 to calculate 
pressure and temperature at specified output times (sec) ) 

Col. 1 Col. 2 Col. 3 

SEP2 (TSPEC(l) ) SEPT2 (TSPEC(l) ) TSPEC(l) 


SEP2 (TSPEC(NF) ) SEPT2 (TSPEC(NF) ) TSPEC(NF) 
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stress. out (through-thickness distributions of stress (Pa) at specified 


output times (sec) ) 


Col. 1 

Col. 2 

Col. 3 

Col. 4 

Col. 5 

Col. 6 

Col. 7 

Col. 8 

Position 


Vl 

&* 

cT 

Vi 

a m 


TSPEC(l) 

(m) 




. 


• 

Position 

/ \ 


cT 

Vs 

cT 

<f 

*2*3 


TSPEC 

(NF) 


(m) 
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temp. out (through-thickness distributions of temperature ( C) at specified 

output times (sec) ) 

Col. 1 Col. 2 Col. 3 

Position (m) T (TSPEC(l) ) TSPEC(l) 


Position (m) T (TSPEC(NF) ) TSPEC(NF) 
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volat.out (through-thickness distributions of moisture content and rate of 
temperature change CC/sec) at specified output times (sec) ) 


Col. 1 

Position (m) 

Col. 2 
MC 

(TSPEC(l) ) 

Col. 3 
dr/dt 

(TSPEC(l) ) 

Col. 4 
TSPEC(l) 

Position (m) 

MC 

(TSPEC(NF) ) 

dT/3t 

(TSPEC(NF) ) 

TSPEC(NF) 


193 


APPENDIX H 


INPUT AND OUTPUT FILES OF RTG 

PROGRAM 


All of the input files of the RTG program have the same formats and 
contents as the input files of the ablative composite plate program except for 
the input file named: "input.in." Also, for the RTG program, there is one 
additional input file named: "rtg.in." The RTG program has the same output 
files as the ablative composite plate program. The formats and contents of 
the two input files (input.in and rtg.in) are shown in the following sections. 

H. 1 Formats and Contents of Input File "input.in" 

input.in 

I. TRATESPEC 

Constant heating rate CC/sec) 

2. SEP2 

Number of nodes using IFDM solution scheme for pressure (Eqn. 4-46) 

3. SEPT2 

Number of nodes using 2nd EFDM solution scheme for temperature 
(Eqn. 4-48) 

4. NF 

Number of specified output time 

5. TSPEC(l) 

First specified output time (sec) 
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TSPEC(NF) 

NFth specified output time (sec) 

6. GN 

Number of char volumes at which the permeability of the porous solid 
are specified 

7. KN 

Number of temperatures at which coefficients of thermal conductivity 
are specified 

8. NET 

Number of temperatures at which mechanical properties are specified 

9. NEC 

Number of char volumes at which mechanical properties are specified 

10. CN 

Number of temperatures at which the specific heats of the porous 
virgin and charred solids are specified 

11. MN 

Number of temperatures at which the gas viscosities are specified 

12. TYPEC 

Type of char reaction. For RTG problem enter 4 

13. TYPEW 

Type of evaporation reaction. For RTG problem enter 4 

14. NRTG 

Types of chemical reactions 

15. PRSP 

Number of pressure at which the reaction constants are specified 
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16 NODE 


Number of nodes in the mesh 

H.2 Formats and Contents of Input File "rtg.in" 
rtg.in 

Pi 

First pressure at which the Arrhenius constants are specified (Pa) 


Pprsp 

PRSPth pressure at which the Arrhenius constants are specified (Pa) 

*,(Pi> 

Activation energy (J/mole) for first type chemical reaction evaluated at 
first specified pressure 

^a, (PpRSp) 

Activation energy (J/mole) for first type chemical reaction evaluated at 
PRSPth specified pressure 


Activation energy (J/mole) for NRTGth type chemical reaction 
evaluated at first specified pressure 


Activation energy (J/mole) for the NRTGth type chemical reaction 
evaluated at PRSPth specified pressure 
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